Fix compilation error in debug mode.
[python.git] / Python / dtoa.c
blobc0559b487272eca8a2792cb0fccd21a49df281b2
1 /****************************************************************
3 * The author of this software is David M. Gay.
5 * Copyright (c) 1991, 2000, 2001 by Lucent Technologies.
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.
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.
18 ***************************************************************/
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.
25 * Please remember to check http://www.netlib.org/fp regularly (and especially
26 * before any Python release) for bugfixes and updates.
28 * The major modifications from Gay's original code are as follows:
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.
36 * 1. We use PyMem_Malloc and PyMem_Free in place of malloc and free.
38 * 2. The public functions strtod, dtoa and freedtoa all now have
39 * a _Py_dg_ prefix.
41 * 3. Instead of assuming that PyMem_Malloc always succeeds, we thread
42 * PyMem_Malloc failures through the code. The functions
44 * Balloc, multadd, s2b, i2b, mult, pow5mult, lshift, diff, d2b
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.
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.
57 * 5. The code has been reformatted to better fit with Python's
58 * C style guide (PEP 7).
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.
64 * 7. _Py_dg_strtod has been modified so that it doesn't accept strings with
65 * leading whitespace.
67 ***************************************************************/
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). */
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.
85 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
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).
92 * Inspired loosely by William D. Clinger's paper "How to Read Floating
93 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
95 * Modifications:
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).
115 /* Linking of Python's #defines to Gay's #defines starts here. */
117 #include "Python.h"
119 #define IEEE_8087
121 /* if PY_NO_SHORT_FLOAT_REPR is defined, then don't even try to compile
122 the following code */
123 #ifndef PY_NO_SHORT_FLOAT_REPR
125 #include "float.h"
127 #define MALLOC PyMem_Malloc
128 #define FREE PyMem_Free
130 /* This code should also work for ARM mixed-endian format on little-endian
131 machines, where doubles have byte order 45670123 (in increasing address
132 order, 0 being the least significant byte). */
133 #ifdef DOUBLE_IS_LITTLE_ENDIAN_IEEE754
134 # define IEEE_8087
135 #endif
136 #if defined(DOUBLE_IS_BIG_ENDIAN_IEEE754) || \
137 defined(DOUBLE_IS_ARM_MIXED_ENDIAN_IEEE754)
138 # define IEEE_MC68k
139 #endif
140 #if defined(IEEE_8087) + defined(IEEE_MC68k) != 1
141 #error "Exactly one of IEEE_8087 or IEEE_MC68k should be defined."
142 #endif
144 /* The code below assumes that the endianness of integers matches the
145 endianness of the two 32-bit words of a double. Check this. */
146 #if defined(WORDS_BIGENDIAN) && (defined(DOUBLE_IS_LITTLE_ENDIAN_IEEE754) || \
147 defined(DOUBLE_IS_ARM_MIXED_ENDIAN_IEEE754))
148 #error "doubles and ints have incompatible endianness"
149 #endif
151 #if !defined(WORDS_BIGENDIAN) && defined(DOUBLE_IS_BIG_ENDIAN_IEEE754)
152 #error "doubles and ints have incompatible endianness"
153 #endif
156 #if defined(HAVE_UINT32_T) && defined(HAVE_INT32_T)
157 typedef PY_UINT32_T ULong;
158 typedef PY_INT32_T Long;
159 #else
160 #error "Failed to find an exact-width 32-bit integer type"
161 #endif
163 #if defined(HAVE_UINT64_T)
164 #define ULLong PY_UINT64_T
165 #else
166 #undef ULLong
167 #endif
169 #undef DEBUG
170 #ifdef Py_DEBUG
171 #define DEBUG
172 #endif
174 /* End Python #define linking */
176 #ifdef DEBUG
177 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
178 #endif
180 #ifndef PRIVATE_MEM
181 #define PRIVATE_MEM 2304
182 #endif
183 #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
184 static double private_mem[PRIVATE_mem], *pmem_next = private_mem;
186 #ifdef __cplusplus
187 extern "C" {
188 #endif
190 typedef union { double d; ULong L[2]; } U;
192 #ifdef IEEE_8087
193 #define word0(x) (x)->L[1]
194 #define word1(x) (x)->L[0]
195 #else
196 #define word0(x) (x)->L[0]
197 #define word1(x) (x)->L[1]
198 #endif
199 #define dval(x) (x)->d
201 #ifndef STRTOD_DIGLIM
202 #define STRTOD_DIGLIM 40
203 #endif
205 #ifdef DIGLIM_DEBUG
206 extern int strtod_diglim;
207 #else
208 #define strtod_diglim STRTOD_DIGLIM
209 #endif
211 /* The following definition of Storeinc is appropriate for MIPS processors.
212 * An alternative that might be better on some machines is
213 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
215 #if defined(IEEE_8087)
216 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
217 ((unsigned short *)a)[0] = (unsigned short)c, a++)
218 #else
219 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
220 ((unsigned short *)a)[1] = (unsigned short)c, a++)
221 #endif
223 /* #define P DBL_MANT_DIG */
224 /* Ten_pmax = floor(P*log(2)/log(5)) */
225 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
226 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
227 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
229 #define Exp_shift 20
230 #define Exp_shift1 20
231 #define Exp_msk1 0x100000
232 #define Exp_msk11 0x100000
233 #define Exp_mask 0x7ff00000
234 #define P 53
235 #define Nbits 53
236 #define Bias 1023
237 #define Emax 1023
238 #define Emin (-1022)
239 #define Exp_1 0x3ff00000
240 #define Exp_11 0x3ff00000
241 #define Ebits 11
242 #define Frac_mask 0xfffff
243 #define Frac_mask1 0xfffff
244 #define Ten_pmax 22
245 #define Bletch 0x10
246 #define Bndry_mask 0xfffff
247 #define Bndry_mask1 0xfffff
248 #define LSB 1
249 #define Sign_bit 0x80000000
250 #define Log2P 1
251 #define Tiny0 0
252 #define Tiny1 1
253 #define Quick_max 14
254 #define Int_max 14
256 #ifndef Flt_Rounds
257 #ifdef FLT_ROUNDS
258 #define Flt_Rounds FLT_ROUNDS
259 #else
260 #define Flt_Rounds 1
261 #endif
262 #endif /*Flt_Rounds*/
264 #define Rounding Flt_Rounds
266 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
267 #define Big1 0xffffffff
269 /* struct BCinfo is used to pass information from _Py_dg_strtod to bigcomp */
271 typedef struct BCinfo BCinfo;
272 struct
273 BCinfo {
274 int dp0, dp1, dplen, dsign, e0, inexact;
275 int nd, nd0, rounding, scale, uflchk;
278 #define FFFFFFFF 0xffffffffUL
280 #define Kmax 7
282 /* struct Bigint is used to represent arbitrary-precision integers. These
283 integers are stored in sign-magnitude format, with the magnitude stored as
284 an array of base 2**32 digits. Bigints are always normalized: if x is a
285 Bigint then x->wds >= 1, and either x->wds == 1 or x[wds-1] is nonzero.
287 The Bigint fields are as follows:
289 - next is a header used by Balloc and Bfree to keep track of lists
290 of freed Bigints; it's also used for the linked list of
291 powers of 5 of the form 5**2**i used by pow5mult.
292 - k indicates which pool this Bigint was allocated from
293 - maxwds is the maximum number of words space was allocated for
294 (usually maxwds == 2**k)
295 - sign is 1 for negative Bigints, 0 for positive. The sign is unused
296 (ignored on inputs, set to 0 on outputs) in almost all operations
297 involving Bigints: a notable exception is the diff function, which
298 ignores signs on inputs but sets the sign of the output correctly.
299 - wds is the actual number of significant words
300 - x contains the vector of words (digits) for this Bigint, from least
301 significant (x[0]) to most significant (x[wds-1]).
304 struct
305 Bigint {
306 struct Bigint *next;
307 int k, maxwds, sign, wds;
308 ULong x[1];
311 typedef struct Bigint Bigint;
313 /* Memory management: memory is allocated from, and returned to, Kmax+1 pools
314 of memory, where pool k (0 <= k <= Kmax) is for Bigints b with b->maxwds ==
315 1 << k. These pools are maintained as linked lists, with freelist[k]
316 pointing to the head of the list for pool k.
318 On allocation, if there's no free slot in the appropriate pool, MALLOC is
319 called to get more memory. This memory is not returned to the system until
320 Python quits. There's also a private memory pool that's allocated from
321 in preference to using MALLOC.
323 For Bigints with more than (1 << Kmax) digits (which implies at least 1233
324 decimal digits), memory is directly allocated using MALLOC, and freed using
325 FREE.
327 XXX: it would be easy to bypass this memory-management system and
328 translate each call to Balloc into a call to PyMem_Malloc, and each
329 Bfree to PyMem_Free. Investigate whether this has any significant
330 performance on impact. */
332 static Bigint *freelist[Kmax+1];
334 /* Allocate space for a Bigint with up to 1<<k digits */
336 static Bigint *
337 Balloc(int k)
339 int x;
340 Bigint *rv;
341 unsigned int len;
343 if (k <= Kmax && (rv = freelist[k]))
344 freelist[k] = rv->next;
345 else {
346 x = 1 << k;
347 len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
348 /sizeof(double);
349 if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) {
350 rv = (Bigint*)pmem_next;
351 pmem_next += len;
353 else {
354 rv = (Bigint*)MALLOC(len*sizeof(double));
355 if (rv == NULL)
356 return NULL;
358 rv->k = k;
359 rv->maxwds = x;
361 rv->sign = rv->wds = 0;
362 return rv;
365 /* Free a Bigint allocated with Balloc */
367 static void
368 Bfree(Bigint *v)
370 if (v) {
371 if (v->k > Kmax)
372 FREE((void*)v);
373 else {
374 v->next = freelist[v->k];
375 freelist[v->k] = v;
380 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
381 y->wds*sizeof(Long) + 2*sizeof(int))
383 /* Multiply a Bigint b by m and add a. Either modifies b in place and returns
384 a pointer to the modified b, or Bfrees b and returns a pointer to a copy.
385 On failure, return NULL. In this case, b will have been already freed. */
387 static Bigint *
388 multadd(Bigint *b, int m, int a) /* multiply by m and add a */
390 int i, wds;
391 #ifdef ULLong
392 ULong *x;
393 ULLong carry, y;
394 #else
395 ULong carry, *x, y;
396 ULong xi, z;
397 #endif
398 Bigint *b1;
400 wds = b->wds;
401 x = b->x;
402 i = 0;
403 carry = a;
404 do {
405 #ifdef ULLong
406 y = *x * (ULLong)m + carry;
407 carry = y >> 32;
408 *x++ = (ULong)(y & FFFFFFFF);
409 #else
410 xi = *x;
411 y = (xi & 0xffff) * m + carry;
412 z = (xi >> 16) * m + (y >> 16);
413 carry = z >> 16;
414 *x++ = (z << 16) + (y & 0xffff);
415 #endif
417 while(++i < wds);
418 if (carry) {
419 if (wds >= b->maxwds) {
420 b1 = Balloc(b->k+1);
421 if (b1 == NULL){
422 Bfree(b);
423 return NULL;
425 Bcopy(b1, b);
426 Bfree(b);
427 b = b1;
429 b->x[wds++] = (ULong)carry;
430 b->wds = wds;
432 return b;
435 /* convert a string s containing nd decimal digits (possibly containing a
436 decimal separator at position nd0, which is ignored) to a Bigint. This
437 function carries on where the parsing code in _Py_dg_strtod leaves off: on
438 entry, y9 contains the result of converting the first 9 digits. Returns
439 NULL on failure. */
441 static Bigint *
442 s2b(const char *s, int nd0, int nd, ULong y9, int dplen)
444 Bigint *b;
445 int i, k;
446 Long x, y;
448 x = (nd + 8) / 9;
449 for(k = 0, y = 1; x > y; y <<= 1, k++) ;
450 b = Balloc(k);
451 if (b == NULL)
452 return NULL;
453 b->x[0] = y9;
454 b->wds = 1;
456 i = 9;
457 if (9 < nd0) {
458 s += 9;
459 do {
460 b = multadd(b, 10, *s++ - '0');
461 if (b == NULL)
462 return NULL;
463 } while(++i < nd0);
464 s += dplen;
466 else
467 s += dplen + 9;
468 for(; i < nd; i++) {
469 b = multadd(b, 10, *s++ - '0');
470 if (b == NULL)
471 return NULL;
473 return b;
476 /* count leading 0 bits in the 32-bit integer x. */
478 static int
479 hi0bits(ULong x)
481 int k = 0;
483 if (!(x & 0xffff0000)) {
484 k = 16;
485 x <<= 16;
487 if (!(x & 0xff000000)) {
488 k += 8;
489 x <<= 8;
491 if (!(x & 0xf0000000)) {
492 k += 4;
493 x <<= 4;
495 if (!(x & 0xc0000000)) {
496 k += 2;
497 x <<= 2;
499 if (!(x & 0x80000000)) {
500 k++;
501 if (!(x & 0x40000000))
502 return 32;
504 return k;
507 /* count trailing 0 bits in the 32-bit integer y, and shift y right by that
508 number of bits. */
510 static int
511 lo0bits(ULong *y)
513 int k;
514 ULong x = *y;
516 if (x & 7) {
517 if (x & 1)
518 return 0;
519 if (x & 2) {
520 *y = x >> 1;
521 return 1;
523 *y = x >> 2;
524 return 2;
526 k = 0;
527 if (!(x & 0xffff)) {
528 k = 16;
529 x >>= 16;
531 if (!(x & 0xff)) {
532 k += 8;
533 x >>= 8;
535 if (!(x & 0xf)) {
536 k += 4;
537 x >>= 4;
539 if (!(x & 0x3)) {
540 k += 2;
541 x >>= 2;
543 if (!(x & 1)) {
544 k++;
545 x >>= 1;
546 if (!x)
547 return 32;
549 *y = x;
550 return k;
553 /* convert a small nonnegative integer to a Bigint */
555 static Bigint *
556 i2b(int i)
558 Bigint *b;
560 b = Balloc(1);
561 if (b == NULL)
562 return NULL;
563 b->x[0] = i;
564 b->wds = 1;
565 return b;
568 /* multiply two Bigints. Returns a new Bigint, or NULL on failure. Ignores
569 the signs of a and b. */
571 static Bigint *
572 mult(Bigint *a, Bigint *b)
574 Bigint *c;
575 int k, wa, wb, wc;
576 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
577 ULong y;
578 #ifdef ULLong
579 ULLong carry, z;
580 #else
581 ULong carry, z;
582 ULong z2;
583 #endif
585 if (a->wds < b->wds) {
586 c = a;
587 a = b;
588 b = c;
590 k = a->k;
591 wa = a->wds;
592 wb = b->wds;
593 wc = wa + wb;
594 if (wc > a->maxwds)
595 k++;
596 c = Balloc(k);
597 if (c == NULL)
598 return NULL;
599 for(x = c->x, xa = x + wc; x < xa; x++)
600 *x = 0;
601 xa = a->x;
602 xae = xa + wa;
603 xb = b->x;
604 xbe = xb + wb;
605 xc0 = c->x;
606 #ifdef ULLong
607 for(; xb < xbe; xc0++) {
608 if ((y = *xb++)) {
609 x = xa;
610 xc = xc0;
611 carry = 0;
612 do {
613 z = *x++ * (ULLong)y + *xc + carry;
614 carry = z >> 32;
615 *xc++ = (ULong)(z & FFFFFFFF);
617 while(x < xae);
618 *xc = (ULong)carry;
621 #else
622 for(; xb < xbe; xb++, xc0++) {
623 if (y = *xb & 0xffff) {
624 x = xa;
625 xc = xc0;
626 carry = 0;
627 do {
628 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
629 carry = z >> 16;
630 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
631 carry = z2 >> 16;
632 Storeinc(xc, z2, z);
634 while(x < xae);
635 *xc = carry;
637 if (y = *xb >> 16) {
638 x = xa;
639 xc = xc0;
640 carry = 0;
641 z2 = *xc;
642 do {
643 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
644 carry = z >> 16;
645 Storeinc(xc, z, z2);
646 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
647 carry = z2 >> 16;
649 while(x < xae);
650 *xc = z2;
653 #endif
654 for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
655 c->wds = wc;
656 return c;
659 /* p5s is a linked list of powers of 5 of the form 5**(2**i), i >= 2 */
661 static Bigint *p5s;
663 /* multiply the Bigint b by 5**k. Returns a pointer to the result, or NULL on
664 failure; if the returned pointer is distinct from b then the original
665 Bigint b will have been Bfree'd. Ignores the sign of b. */
667 static Bigint *
668 pow5mult(Bigint *b, int k)
670 Bigint *b1, *p5, *p51;
671 int i;
672 static int p05[3] = { 5, 25, 125 };
674 if ((i = k & 3)) {
675 b = multadd(b, p05[i-1], 0);
676 if (b == NULL)
677 return NULL;
680 if (!(k >>= 2))
681 return b;
682 p5 = p5s;
683 if (!p5) {
684 /* first time */
685 p5 = i2b(625);
686 if (p5 == NULL) {
687 Bfree(b);
688 return NULL;
690 p5s = p5;
691 p5->next = 0;
693 for(;;) {
694 if (k & 1) {
695 b1 = mult(b, p5);
696 Bfree(b);
697 b = b1;
698 if (b == NULL)
699 return NULL;
701 if (!(k >>= 1))
702 break;
703 p51 = p5->next;
704 if (!p51) {
705 p51 = mult(p5,p5);
706 if (p51 == NULL) {
707 Bfree(b);
708 return NULL;
710 p51->next = 0;
711 p5->next = p51;
713 p5 = p51;
715 return b;
718 /* shift a Bigint b left by k bits. Return a pointer to the shifted result,
719 or NULL on failure. If the returned pointer is distinct from b then the
720 original b will have been Bfree'd. Ignores the sign of b. */
722 static Bigint *
723 lshift(Bigint *b, int k)
725 int i, k1, n, n1;
726 Bigint *b1;
727 ULong *x, *x1, *xe, z;
729 n = k >> 5;
730 k1 = b->k;
731 n1 = n + b->wds + 1;
732 for(i = b->maxwds; n1 > i; i <<= 1)
733 k1++;
734 b1 = Balloc(k1);
735 if (b1 == NULL) {
736 Bfree(b);
737 return NULL;
739 x1 = b1->x;
740 for(i = 0; i < n; i++)
741 *x1++ = 0;
742 x = b->x;
743 xe = x + b->wds;
744 if (k &= 0x1f) {
745 k1 = 32 - k;
746 z = 0;
747 do {
748 *x1++ = *x << k | z;
749 z = *x++ >> k1;
751 while(x < xe);
752 if ((*x1 = z))
753 ++n1;
755 else do
756 *x1++ = *x++;
757 while(x < xe);
758 b1->wds = n1 - 1;
759 Bfree(b);
760 return b1;
763 /* Do a three-way compare of a and b, returning -1 if a < b, 0 if a == b and
764 1 if a > b. Ignores signs of a and b. */
766 static int
767 cmp(Bigint *a, Bigint *b)
769 ULong *xa, *xa0, *xb, *xb0;
770 int i, j;
772 i = a->wds;
773 j = b->wds;
774 #ifdef DEBUG
775 if (i > 1 && !a->x[i-1])
776 Bug("cmp called with a->x[a->wds-1] == 0");
777 if (j > 1 && !b->x[j-1])
778 Bug("cmp called with b->x[b->wds-1] == 0");
779 #endif
780 if (i -= j)
781 return i;
782 xa0 = a->x;
783 xa = xa0 + j;
784 xb0 = b->x;
785 xb = xb0 + j;
786 for(;;) {
787 if (*--xa != *--xb)
788 return *xa < *xb ? -1 : 1;
789 if (xa <= xa0)
790 break;
792 return 0;
795 /* Take the difference of Bigints a and b, returning a new Bigint. Returns
796 NULL on failure. The signs of a and b are ignored, but the sign of the
797 result is set appropriately. */
799 static Bigint *
800 diff(Bigint *a, Bigint *b)
802 Bigint *c;
803 int i, wa, wb;
804 ULong *xa, *xae, *xb, *xbe, *xc;
805 #ifdef ULLong
806 ULLong borrow, y;
807 #else
808 ULong borrow, y;
809 ULong z;
810 #endif
812 i = cmp(a,b);
813 if (!i) {
814 c = Balloc(0);
815 if (c == NULL)
816 return NULL;
817 c->wds = 1;
818 c->x[0] = 0;
819 return c;
821 if (i < 0) {
822 c = a;
823 a = b;
824 b = c;
825 i = 1;
827 else
828 i = 0;
829 c = Balloc(a->k);
830 if (c == NULL)
831 return NULL;
832 c->sign = i;
833 wa = a->wds;
834 xa = a->x;
835 xae = xa + wa;
836 wb = b->wds;
837 xb = b->x;
838 xbe = xb + wb;
839 xc = c->x;
840 borrow = 0;
841 #ifdef ULLong
842 do {
843 y = (ULLong)*xa++ - *xb++ - borrow;
844 borrow = y >> 32 & (ULong)1;
845 *xc++ = (ULong)(y & FFFFFFFF);
847 while(xb < xbe);
848 while(xa < xae) {
849 y = *xa++ - borrow;
850 borrow = y >> 32 & (ULong)1;
851 *xc++ = (ULong)(y & FFFFFFFF);
853 #else
854 do {
855 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
856 borrow = (y & 0x10000) >> 16;
857 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
858 borrow = (z & 0x10000) >> 16;
859 Storeinc(xc, z, y);
861 while(xb < xbe);
862 while(xa < xae) {
863 y = (*xa & 0xffff) - borrow;
864 borrow = (y & 0x10000) >> 16;
865 z = (*xa++ >> 16) - borrow;
866 borrow = (z & 0x10000) >> 16;
867 Storeinc(xc, z, y);
869 #endif
870 while(!*--xc)
871 wa--;
872 c->wds = wa;
873 return c;
876 /* Given a positive normal double x, return the difference between x and the next
877 double up. Doesn't give correct results for subnormals. */
879 static double
880 ulp(U *x)
882 Long L;
883 U u;
885 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
886 word0(&u) = L;
887 word1(&u) = 0;
888 return dval(&u);
891 /* Convert a Bigint to a double plus an exponent */
893 static double
894 b2d(Bigint *a, int *e)
896 ULong *xa, *xa0, w, y, z;
897 int k;
898 U d;
900 xa0 = a->x;
901 xa = xa0 + a->wds;
902 y = *--xa;
903 #ifdef DEBUG
904 if (!y) Bug("zero y in b2d");
905 #endif
906 k = hi0bits(y);
907 *e = 32 - k;
908 if (k < Ebits) {
909 word0(&d) = Exp_1 | y >> (Ebits - k);
910 w = xa > xa0 ? *--xa : 0;
911 word1(&d) = y << ((32-Ebits) + k) | w >> (Ebits - k);
912 goto ret_d;
914 z = xa > xa0 ? *--xa : 0;
915 if (k -= Ebits) {
916 word0(&d) = Exp_1 | y << k | z >> (32 - k);
917 y = xa > xa0 ? *--xa : 0;
918 word1(&d) = z << k | y >> (32 - k);
920 else {
921 word0(&d) = Exp_1 | y;
922 word1(&d) = z;
924 ret_d:
925 return dval(&d);
928 /* Convert a double to a Bigint plus an exponent. Return NULL on failure.
930 Given a finite nonzero double d, return an odd Bigint b and exponent *e
931 such that fabs(d) = b * 2**e. On return, *bbits gives the number of
932 significant bits of e; that is, 2**(*bbits-1) <= b < 2**(*bbits).
934 If d is zero, then b == 0, *e == -1010, *bbits = 0.
938 static Bigint *
939 d2b(U *d, int *e, int *bits)
941 Bigint *b;
942 int de, k;
943 ULong *x, y, z;
944 int i;
946 b = Balloc(1);
947 if (b == NULL)
948 return NULL;
949 x = b->x;
951 z = word0(d) & Frac_mask;
952 word0(d) &= 0x7fffffff; /* clear sign bit, which we ignore */
953 if ((de = (int)(word0(d) >> Exp_shift)))
954 z |= Exp_msk1;
955 if ((y = word1(d))) {
956 if ((k = lo0bits(&y))) {
957 x[0] = y | z << (32 - k);
958 z >>= k;
960 else
961 x[0] = y;
963 b->wds = (x[1] = z) ? 2 : 1;
965 else {
966 k = lo0bits(&z);
967 x[0] = z;
969 b->wds = 1;
970 k += 32;
972 if (de) {
973 *e = de - Bias - (P-1) + k;
974 *bits = P - k;
976 else {
977 *e = de - Bias - (P-1) + 1 + k;
978 *bits = 32*i - hi0bits(x[i-1]);
980 return b;
983 /* Compute the ratio of two Bigints, as a double. The result may have an
984 error of up to 2.5 ulps. */
986 static double
987 ratio(Bigint *a, Bigint *b)
989 U da, db;
990 int k, ka, kb;
992 dval(&da) = b2d(a, &ka);
993 dval(&db) = b2d(b, &kb);
994 k = ka - kb + 32*(a->wds - b->wds);
995 if (k > 0)
996 word0(&da) += k*Exp_msk1;
997 else {
998 k = -k;
999 word0(&db) += k*Exp_msk1;
1001 return dval(&da) / dval(&db);
1004 static const double
1005 tens[] = {
1006 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1007 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1008 1e20, 1e21, 1e22
1011 static const double
1012 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1013 static const double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
1014 9007199254740992.*9007199254740992.e-256
1015 /* = 2^106 * 1e-256 */
1017 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
1018 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */
1019 #define Scale_Bit 0x10
1020 #define n_bigtens 5
1022 #define ULbits 32
1023 #define kshift 5
1024 #define kmask 31
1027 static int
1028 dshift(Bigint *b, int p2)
1030 int rv = hi0bits(b->x[b->wds-1]) - 4;
1031 if (p2 > 0)
1032 rv -= p2;
1033 return rv & kmask;
1036 /* special case of Bigint division. The quotient is always in the range 0 <=
1037 quotient < 10, and on entry the divisor S is normalized so that its top 4
1038 bits (28--31) are zero and bit 27 is set. */
1040 static int
1041 quorem(Bigint *b, Bigint *S)
1043 int n;
1044 ULong *bx, *bxe, q, *sx, *sxe;
1045 #ifdef ULLong
1046 ULLong borrow, carry, y, ys;
1047 #else
1048 ULong borrow, carry, y, ys;
1049 ULong si, z, zs;
1050 #endif
1052 n = S->wds;
1053 #ifdef DEBUG
1054 /*debug*/ if (b->wds > n)
1055 /*debug*/ Bug("oversize b in quorem");
1056 #endif
1057 if (b->wds < n)
1058 return 0;
1059 sx = S->x;
1060 sxe = sx + --n;
1061 bx = b->x;
1062 bxe = bx + n;
1063 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
1064 #ifdef DEBUG
1065 /*debug*/ if (q > 9)
1066 /*debug*/ Bug("oversized quotient in quorem");
1067 #endif
1068 if (q) {
1069 borrow = 0;
1070 carry = 0;
1071 do {
1072 #ifdef ULLong
1073 ys = *sx++ * (ULLong)q + carry;
1074 carry = ys >> 32;
1075 y = *bx - (ys & FFFFFFFF) - borrow;
1076 borrow = y >> 32 & (ULong)1;
1077 *bx++ = (ULong)(y & FFFFFFFF);
1078 #else
1079 si = *sx++;
1080 ys = (si & 0xffff) * q + carry;
1081 zs = (si >> 16) * q + (ys >> 16);
1082 carry = zs >> 16;
1083 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
1084 borrow = (y & 0x10000) >> 16;
1085 z = (*bx >> 16) - (zs & 0xffff) - borrow;
1086 borrow = (z & 0x10000) >> 16;
1087 Storeinc(bx, z, y);
1088 #endif
1090 while(sx <= sxe);
1091 if (!*bxe) {
1092 bx = b->x;
1093 while(--bxe > bx && !*bxe)
1094 --n;
1095 b->wds = n;
1098 if (cmp(b, S) >= 0) {
1099 q++;
1100 borrow = 0;
1101 carry = 0;
1102 bx = b->x;
1103 sx = S->x;
1104 do {
1105 #ifdef ULLong
1106 ys = *sx++ + carry;
1107 carry = ys >> 32;
1108 y = *bx - (ys & FFFFFFFF) - borrow;
1109 borrow = y >> 32 & (ULong)1;
1110 *bx++ = (ULong)(y & FFFFFFFF);
1111 #else
1112 si = *sx++;
1113 ys = (si & 0xffff) + carry;
1114 zs = (si >> 16) + (ys >> 16);
1115 carry = zs >> 16;
1116 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
1117 borrow = (y & 0x10000) >> 16;
1118 z = (*bx >> 16) - (zs & 0xffff) - borrow;
1119 borrow = (z & 0x10000) >> 16;
1120 Storeinc(bx, z, y);
1121 #endif
1123 while(sx <= sxe);
1124 bx = b->x;
1125 bxe = bx + n;
1126 if (!*bxe) {
1127 while(--bxe > bx && !*bxe)
1128 --n;
1129 b->wds = n;
1132 return q;
1136 /* return 0 on success, -1 on failure */
1138 static int
1139 bigcomp(U *rv, const char *s0, BCinfo *bc)
1141 Bigint *b, *d;
1142 int b2, bbits, d2, dd, dig, dsign, i, j, nd, nd0, p2, p5, speccase;
1144 dsign = bc->dsign;
1145 nd = bc->nd;
1146 nd0 = bc->nd0;
1147 p5 = nd + bc->e0 - 1;
1148 speccase = 0;
1149 if (rv->d == 0.) { /* special case: value near underflow-to-zero */
1150 /* threshold was rounded to zero */
1151 b = i2b(1);
1152 if (b == NULL)
1153 return -1;
1154 p2 = Emin - P + 1;
1155 bbits = 1;
1156 word0(rv) = (P+2) << Exp_shift;
1157 i = 0;
1159 speccase = 1;
1160 --p2;
1161 dsign = 0;
1162 goto have_i;
1165 else
1167 b = d2b(rv, &p2, &bbits);
1168 if (b == NULL)
1169 return -1;
1171 p2 -= bc->scale;
1172 /* floor(log2(rv)) == bbits - 1 + p2 */
1173 /* Check for denormal case. */
1174 i = P - bbits;
1175 if (i > (j = P - Emin - 1 + p2)) {
1176 i = j;
1179 b = lshift(b, ++i);
1180 if (b == NULL)
1181 return -1;
1182 b->x[0] |= 1;
1184 have_i:
1185 p2 -= p5 + i;
1186 d = i2b(1);
1187 if (d == NULL) {
1188 Bfree(b);
1189 return -1;
1191 /* Arrange for convenient computation of quotients:
1192 * shift left if necessary so divisor has 4 leading 0 bits.
1194 if (p5 > 0) {
1195 d = pow5mult(d, p5);
1196 if (d == NULL) {
1197 Bfree(b);
1198 return -1;
1201 else if (p5 < 0) {
1202 b = pow5mult(b, -p5);
1203 if (b == NULL) {
1204 Bfree(d);
1205 return -1;
1208 if (p2 > 0) {
1209 b2 = p2;
1210 d2 = 0;
1212 else {
1213 b2 = 0;
1214 d2 = -p2;
1216 i = dshift(d, d2);
1217 if ((b2 += i) > 0) {
1218 b = lshift(b, b2);
1219 if (b == NULL) {
1220 Bfree(d);
1221 return -1;
1224 if ((d2 += i) > 0) {
1225 d = lshift(d, d2);
1226 if (d == NULL) {
1227 Bfree(b);
1228 return -1;
1232 /* Now b/d = exactly half-way between the two floating-point values */
1233 /* on either side of the input string. Compute first digit of b/d. */
1235 if (!(dig = quorem(b,d))) {
1236 b = multadd(b, 10, 0); /* very unlikely */
1237 if (b == NULL) {
1238 Bfree(d);
1239 return -1;
1241 dig = quorem(b,d);
1244 /* Compare b/d with s0 */
1246 assert(nd > 0);
1247 dd = 9999; /* silence gcc compiler warning */
1248 for(i = 0; i < nd0; ) {
1249 if ((dd = s0[i++] - '0' - dig))
1250 goto ret;
1251 if (!b->x[0] && b->wds == 1) {
1252 if (i < nd)
1253 dd = 1;
1254 goto ret;
1256 b = multadd(b, 10, 0);
1257 if (b == NULL) {
1258 Bfree(d);
1259 return -1;
1261 dig = quorem(b,d);
1263 for(j = bc->dp1; i++ < nd;) {
1264 if ((dd = s0[j++] - '0' - dig))
1265 goto ret;
1266 if (!b->x[0] && b->wds == 1) {
1267 if (i < nd)
1268 dd = 1;
1269 goto ret;
1271 b = multadd(b, 10, 0);
1272 if (b == NULL) {
1273 Bfree(d);
1274 return -1;
1276 dig = quorem(b,d);
1278 if (b->x[0] || b->wds > 1)
1279 dd = -1;
1280 ret:
1281 Bfree(b);
1282 Bfree(d);
1283 if (speccase) {
1284 if (dd <= 0)
1285 rv->d = 0.;
1287 else if (dd < 0) {
1288 if (!dsign) /* does not happen for round-near */
1289 retlow1:
1290 dval(rv) -= ulp(rv);
1292 else if (dd > 0) {
1293 if (dsign) {
1294 rethi1:
1295 dval(rv) += ulp(rv);
1298 else {
1299 /* Exact half-way case: apply round-even rule. */
1300 if (word1(rv) & 1) {
1301 if (dsign)
1302 goto rethi1;
1303 goto retlow1;
1307 return 0;
1310 double
1311 _Py_dg_strtod(const char *s00, char **se)
1313 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, e, e1, error;
1314 int esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
1315 const char *s, *s0, *s1;
1316 double aadj, aadj1;
1317 Long L;
1318 U aadj2, adj, rv, rv0;
1319 ULong y, z;
1320 BCinfo bc;
1321 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
1323 sign = nz0 = nz = bc.dplen = bc.uflchk = 0;
1324 dval(&rv) = 0.;
1325 for(s = s00;;s++) switch(*s) {
1326 case '-':
1327 sign = 1;
1328 /* no break */
1329 case '+':
1330 if (*++s)
1331 goto break2;
1332 /* no break */
1333 case 0:
1334 goto ret0;
1335 /* modify original dtoa.c so that it doesn't accept leading whitespace
1336 case '\t':
1337 case '\n':
1338 case '\v':
1339 case '\f':
1340 case '\r':
1341 case ' ':
1342 continue;
1344 default:
1345 goto break2;
1347 break2:
1348 if (*s == '0') {
1349 nz0 = 1;
1350 while(*++s == '0') ;
1351 if (!*s)
1352 goto ret;
1354 s0 = s;
1355 y = z = 0;
1356 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
1357 if (nd < 9)
1358 y = 10*y + c - '0';
1359 else if (nd < 16)
1360 z = 10*z + c - '0';
1361 nd0 = nd;
1362 bc.dp0 = bc.dp1 = s - s0;
1363 if (c == '.') {
1364 c = *++s;
1365 bc.dp1 = s - s0;
1366 bc.dplen = bc.dp1 - bc.dp0;
1367 if (!nd) {
1368 for(; c == '0'; c = *++s)
1369 nz++;
1370 if (c > '0' && c <= '9') {
1371 s0 = s;
1372 nf += nz;
1373 nz = 0;
1374 goto have_dig;
1376 goto dig_done;
1378 for(; c >= '0' && c <= '9'; c = *++s) {
1379 have_dig:
1380 nz++;
1381 if (c -= '0') {
1382 nf += nz;
1383 for(i = 1; i < nz; i++)
1384 if (nd++ < 9)
1385 y *= 10;
1386 else if (nd <= DBL_DIG + 1)
1387 z *= 10;
1388 if (nd++ < 9)
1389 y = 10*y + c;
1390 else if (nd <= DBL_DIG + 1)
1391 z = 10*z + c;
1392 nz = 0;
1396 dig_done:
1397 e = 0;
1398 if (c == 'e' || c == 'E') {
1399 if (!nd && !nz && !nz0) {
1400 goto ret0;
1402 s00 = s;
1403 esign = 0;
1404 switch(c = *++s) {
1405 case '-':
1406 esign = 1;
1407 case '+':
1408 c = *++s;
1410 if (c >= '0' && c <= '9') {
1411 while(c == '0')
1412 c = *++s;
1413 if (c > '0' && c <= '9') {
1414 L = c - '0';
1415 s1 = s;
1416 while((c = *++s) >= '0' && c <= '9')
1417 L = 10*L + c - '0';
1418 if (s - s1 > 8 || L > 19999)
1419 /* Avoid confusion from exponents
1420 * so large that e might overflow.
1422 e = 19999; /* safe for 16 bit ints */
1423 else
1424 e = (int)L;
1425 if (esign)
1426 e = -e;
1428 else
1429 e = 0;
1431 else
1432 s = s00;
1434 if (!nd) {
1435 if (!nz && !nz0) {
1436 ret0:
1437 s = s00;
1438 sign = 0;
1440 goto ret;
1442 bc.e0 = e1 = e -= nf;
1444 /* Now we have nd0 digits, starting at s0, followed by a
1445 * decimal point, followed by nd-nd0 digits. The number we're
1446 * after is the integer represented by those digits times
1447 * 10**e */
1449 if (!nd0)
1450 nd0 = nd;
1451 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
1452 dval(&rv) = y;
1453 if (k > 9) {
1454 dval(&rv) = tens[k - 9] * dval(&rv) + z;
1456 bd0 = 0;
1457 if (nd <= DBL_DIG
1458 && Flt_Rounds == 1
1460 if (!e)
1461 goto ret;
1462 if (e > 0) {
1463 if (e <= Ten_pmax) {
1464 dval(&rv) *= tens[e];
1465 goto ret;
1467 i = DBL_DIG - nd;
1468 if (e <= Ten_pmax + i) {
1469 /* A fancier test would sometimes let us do
1470 * this for larger i values.
1472 e -= i;
1473 dval(&rv) *= tens[i];
1474 dval(&rv) *= tens[e];
1475 goto ret;
1478 else if (e >= -Ten_pmax) {
1479 dval(&rv) /= tens[-e];
1480 goto ret;
1483 e1 += nd - k;
1485 bc.scale = 0;
1487 /* Get starting approximation = rv * 10**e1 */
1489 if (e1 > 0) {
1490 if ((i = e1 & 15))
1491 dval(&rv) *= tens[i];
1492 if (e1 &= ~15) {
1493 if (e1 > DBL_MAX_10_EXP) {
1494 ovfl:
1495 errno = ERANGE;
1496 /* Can't trust HUGE_VAL */
1497 word0(&rv) = Exp_mask;
1498 word1(&rv) = 0;
1499 goto ret;
1501 e1 >>= 4;
1502 for(j = 0; e1 > 1; j++, e1 >>= 1)
1503 if (e1 & 1)
1504 dval(&rv) *= bigtens[j];
1505 /* The last multiplication could overflow. */
1506 word0(&rv) -= P*Exp_msk1;
1507 dval(&rv) *= bigtens[j];
1508 if ((z = word0(&rv) & Exp_mask)
1509 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
1510 goto ovfl;
1511 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
1512 /* set to largest number */
1513 /* (Can't trust DBL_MAX) */
1514 word0(&rv) = Big0;
1515 word1(&rv) = Big1;
1517 else
1518 word0(&rv) += P*Exp_msk1;
1521 else if (e1 < 0) {
1522 e1 = -e1;
1523 if ((i = e1 & 15))
1524 dval(&rv) /= tens[i];
1525 if (e1 >>= 4) {
1526 if (e1 >= 1 << n_bigtens)
1527 goto undfl;
1528 if (e1 & Scale_Bit)
1529 bc.scale = 2*P;
1530 for(j = 0; e1 > 0; j++, e1 >>= 1)
1531 if (e1 & 1)
1532 dval(&rv) *= tinytens[j];
1533 if (bc.scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask)
1534 >> Exp_shift)) > 0) {
1535 /* scaled rv is denormal; clear j low bits */
1536 if (j >= 32) {
1537 word1(&rv) = 0;
1538 if (j >= 53)
1539 word0(&rv) = (P+2)*Exp_msk1;
1540 else
1541 word0(&rv) &= 0xffffffff << (j-32);
1543 else
1544 word1(&rv) &= 0xffffffff << j;
1546 if (!dval(&rv)) {
1547 undfl:
1548 dval(&rv) = 0.;
1549 errno = ERANGE;
1550 goto ret;
1555 /* Now the hard part -- adjusting rv to the correct value.*/
1557 /* Put digits into bd: true value = bd * 10^e */
1559 bc.nd = nd;
1560 bc.nd0 = nd0; /* Only needed if nd > strtod_diglim, but done here */
1561 /* to silence an erroneous warning about bc.nd0 */
1562 /* possibly not being initialized. */
1563 if (nd > strtod_diglim) {
1564 /* ASSERT(strtod_diglim >= 18); 18 == one more than the */
1565 /* minimum number of decimal digits to distinguish double values */
1566 /* in IEEE arithmetic. */
1567 i = j = 18;
1568 if (i > nd0)
1569 j += bc.dplen;
1570 for(;;) {
1571 if (--j <= bc.dp1 && j >= bc.dp0)
1572 j = bc.dp0 - 1;
1573 if (s0[j] != '0')
1574 break;
1575 --i;
1577 e += nd - i;
1578 nd = i;
1579 if (nd0 > nd)
1580 nd0 = nd;
1581 if (nd < 9) { /* must recompute y */
1582 y = 0;
1583 for(i = 0; i < nd0; ++i)
1584 y = 10*y + s0[i] - '0';
1585 for(j = bc.dp1; i < nd; ++i)
1586 y = 10*y + s0[j++] - '0';
1589 bd0 = s2b(s0, nd0, nd, y, bc.dplen);
1590 if (bd0 == NULL)
1591 goto failed_malloc;
1593 for(;;) {
1594 bd = Balloc(bd0->k);
1595 if (bd == NULL) {
1596 Bfree(bd0);
1597 goto failed_malloc;
1599 Bcopy(bd, bd0);
1600 bb = d2b(&rv, &bbe, &bbbits); /* rv = bb * 2^bbe */
1601 if (bb == NULL) {
1602 Bfree(bd);
1603 Bfree(bd0);
1604 goto failed_malloc;
1606 bs = i2b(1);
1607 if (bs == NULL) {
1608 Bfree(bb);
1609 Bfree(bd);
1610 Bfree(bd0);
1611 goto failed_malloc;
1614 if (e >= 0) {
1615 bb2 = bb5 = 0;
1616 bd2 = bd5 = e;
1618 else {
1619 bb2 = bb5 = -e;
1620 bd2 = bd5 = 0;
1622 if (bbe >= 0)
1623 bb2 += bbe;
1624 else
1625 bd2 -= bbe;
1626 bs2 = bb2;
1627 j = bbe - bc.scale;
1628 i = j + bbbits - 1; /* logb(rv) */
1629 if (i < Emin) /* denormal */
1630 j += P - Emin;
1631 else
1632 j = P + 1 - bbbits;
1633 bb2 += j;
1634 bd2 += j;
1635 bd2 += bc.scale;
1636 i = bb2 < bd2 ? bb2 : bd2;
1637 if (i > bs2)
1638 i = bs2;
1639 if (i > 0) {
1640 bb2 -= i;
1641 bd2 -= i;
1642 bs2 -= i;
1644 if (bb5 > 0) {
1645 bs = pow5mult(bs, bb5);
1646 if (bs == NULL) {
1647 Bfree(bb);
1648 Bfree(bd);
1649 Bfree(bd0);
1650 goto failed_malloc;
1652 bb1 = mult(bs, bb);
1653 Bfree(bb);
1654 bb = bb1;
1655 if (bb == NULL) {
1656 Bfree(bs);
1657 Bfree(bd);
1658 Bfree(bd0);
1659 goto failed_malloc;
1662 if (bb2 > 0) {
1663 bb = lshift(bb, bb2);
1664 if (bb == NULL) {
1665 Bfree(bs);
1666 Bfree(bd);
1667 Bfree(bd0);
1668 goto failed_malloc;
1671 if (bd5 > 0) {
1672 bd = pow5mult(bd, bd5);
1673 if (bd == NULL) {
1674 Bfree(bb);
1675 Bfree(bs);
1676 Bfree(bd0);
1677 goto failed_malloc;
1680 if (bd2 > 0) {
1681 bd = lshift(bd, bd2);
1682 if (bd == NULL) {
1683 Bfree(bb);
1684 Bfree(bs);
1685 Bfree(bd0);
1686 goto failed_malloc;
1689 if (bs2 > 0) {
1690 bs = lshift(bs, bs2);
1691 if (bs == NULL) {
1692 Bfree(bb);
1693 Bfree(bd);
1694 Bfree(bd0);
1695 goto failed_malloc;
1698 delta = diff(bb, bd);
1699 if (delta == NULL) {
1700 Bfree(bb);
1701 Bfree(bs);
1702 Bfree(bd);
1703 Bfree(bd0);
1704 goto failed_malloc;
1706 bc.dsign = delta->sign;
1707 delta->sign = 0;
1708 i = cmp(delta, bs);
1709 if (bc.nd > nd && i <= 0) {
1710 if (bc.dsign)
1711 break; /* Must use bigcomp(). */
1713 bc.nd = nd;
1714 i = -1; /* Discarded digits make delta smaller. */
1718 if (i < 0) {
1719 /* Error is less than half an ulp -- check for
1720 * special case of mantissa a power of two.
1722 if (bc.dsign || word1(&rv) || word0(&rv) & Bndry_mask
1723 || (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1
1725 break;
1727 if (!delta->x[0] && delta->wds <= 1) {
1728 /* exact result */
1729 break;
1731 delta = lshift(delta,Log2P);
1732 if (delta == NULL) {
1733 Bfree(bb);
1734 Bfree(bs);
1735 Bfree(bd);
1736 Bfree(bd0);
1737 goto failed_malloc;
1739 if (cmp(delta, bs) > 0)
1740 goto drop_down;
1741 break;
1743 if (i == 0) {
1744 /* exactly half-way between */
1745 if (bc.dsign) {
1746 if ((word0(&rv) & Bndry_mask1) == Bndry_mask1
1747 && word1(&rv) == (
1748 (bc.scale &&
1749 (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1) ?
1750 (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
1751 0xffffffff)) {
1752 /*boundary case -- increment exponent*/
1753 word0(&rv) = (word0(&rv) & Exp_mask)
1754 + Exp_msk1
1756 word1(&rv) = 0;
1757 bc.dsign = 0;
1758 break;
1761 else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) {
1762 drop_down:
1763 /* boundary case -- decrement exponent */
1764 if (bc.scale) {
1765 L = word0(&rv) & Exp_mask;
1766 if (L <= (2*P+1)*Exp_msk1) {
1767 if (L > (P+2)*Exp_msk1)
1768 /* round even ==> */
1769 /* accept rv */
1770 break;
1771 /* rv = smallest denormal */
1772 if (bc.nd >nd) {
1773 bc.uflchk = 1;
1774 break;
1776 goto undfl;
1779 L = (word0(&rv) & Exp_mask) - Exp_msk1;
1780 word0(&rv) = L | Bndry_mask1;
1781 word1(&rv) = 0xffffffff;
1782 break;
1784 if (!(word1(&rv) & LSB))
1785 break;
1786 if (bc.dsign)
1787 dval(&rv) += ulp(&rv);
1788 else {
1789 dval(&rv) -= ulp(&rv);
1790 if (!dval(&rv)) {
1791 if (bc.nd >nd) {
1792 bc.uflchk = 1;
1793 break;
1795 goto undfl;
1798 bc.dsign = 1 - bc.dsign;
1799 break;
1801 if ((aadj = ratio(delta, bs)) <= 2.) {
1802 if (bc.dsign)
1803 aadj = aadj1 = 1.;
1804 else if (word1(&rv) || word0(&rv) & Bndry_mask) {
1805 if (word1(&rv) == Tiny1 && !word0(&rv)) {
1806 if (bc.nd >nd) {
1807 bc.uflchk = 1;
1808 break;
1810 goto undfl;
1812 aadj = 1.;
1813 aadj1 = -1.;
1815 else {
1816 /* special case -- power of FLT_RADIX to be */
1817 /* rounded down... */
1819 if (aadj < 2./FLT_RADIX)
1820 aadj = 1./FLT_RADIX;
1821 else
1822 aadj *= 0.5;
1823 aadj1 = -aadj;
1826 else {
1827 aadj *= 0.5;
1828 aadj1 = bc.dsign ? aadj : -aadj;
1829 if (Flt_Rounds == 0)
1830 aadj1 += 0.5;
1832 y = word0(&rv) & Exp_mask;
1834 /* Check for overflow */
1836 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
1837 dval(&rv0) = dval(&rv);
1838 word0(&rv) -= P*Exp_msk1;
1839 adj.d = aadj1 * ulp(&rv);
1840 dval(&rv) += adj.d;
1841 if ((word0(&rv) & Exp_mask) >=
1842 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
1843 if (word0(&rv0) == Big0 && word1(&rv0) == Big1)
1844 goto ovfl;
1845 word0(&rv) = Big0;
1846 word1(&rv) = Big1;
1847 goto cont;
1849 else
1850 word0(&rv) += P*Exp_msk1;
1852 else {
1853 if (bc.scale && y <= 2*P*Exp_msk1) {
1854 if (aadj <= 0x7fffffff) {
1855 if ((z = (ULong)aadj) <= 0)
1856 z = 1;
1857 aadj = z;
1858 aadj1 = bc.dsign ? aadj : -aadj;
1860 dval(&aadj2) = aadj1;
1861 word0(&aadj2) += (2*P+1)*Exp_msk1 - y;
1862 aadj1 = dval(&aadj2);
1864 adj.d = aadj1 * ulp(&rv);
1865 dval(&rv) += adj.d;
1867 z = word0(&rv) & Exp_mask;
1868 if (bc.nd == nd) {
1869 if (!bc.scale)
1870 if (y == z) {
1871 /* Can we stop now? */
1872 L = (Long)aadj;
1873 aadj -= L;
1874 /* The tolerances below are conservative. */
1875 if (bc.dsign || word1(&rv) || word0(&rv) & Bndry_mask) {
1876 if (aadj < .4999999 || aadj > .5000001)
1877 break;
1879 else if (aadj < .4999999/FLT_RADIX)
1880 break;
1883 cont:
1884 Bfree(bb);
1885 Bfree(bd);
1886 Bfree(bs);
1887 Bfree(delta);
1889 Bfree(bb);
1890 Bfree(bd);
1891 Bfree(bs);
1892 Bfree(bd0);
1893 Bfree(delta);
1894 if (bc.nd > nd) {
1895 error = bigcomp(&rv, s0, &bc);
1896 if (error)
1897 goto failed_malloc;
1900 if (bc.scale) {
1901 word0(&rv0) = Exp_1 - 2*P*Exp_msk1;
1902 word1(&rv0) = 0;
1903 dval(&rv) *= dval(&rv0);
1904 /* try to avoid the bug of testing an 8087 register value */
1905 if (!(word0(&rv) & Exp_mask))
1906 errno = ERANGE;
1908 ret:
1909 if (se)
1910 *se = (char *)s;
1911 return sign ? -dval(&rv) : dval(&rv);
1913 failed_malloc:
1914 if (se)
1915 *se = (char *)s00;
1916 errno = ENOMEM;
1917 return -1.0;
1920 static char *
1921 rv_alloc(int i)
1923 int j, k, *r;
1925 j = sizeof(ULong);
1926 for(k = 0;
1927 sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= (unsigned)i;
1928 j <<= 1)
1929 k++;
1930 r = (int*)Balloc(k);
1931 if (r == NULL)
1932 return NULL;
1933 *r = k;
1934 return (char *)(r+1);
1937 static char *
1938 nrv_alloc(char *s, char **rve, int n)
1940 char *rv, *t;
1942 rv = rv_alloc(n);
1943 if (rv == NULL)
1944 return NULL;
1945 t = rv;
1946 while((*t = *s++)) t++;
1947 if (rve)
1948 *rve = t;
1949 return rv;
1952 /* freedtoa(s) must be used to free values s returned by dtoa
1953 * when MULTIPLE_THREADS is #defined. It should be used in all cases,
1954 * but for consistency with earlier versions of dtoa, it is optional
1955 * when MULTIPLE_THREADS is not defined.
1958 void
1959 _Py_dg_freedtoa(char *s)
1961 Bigint *b = (Bigint *)((int *)s - 1);
1962 b->maxwds = 1 << (b->k = *(int*)b);
1963 Bfree(b);
1966 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
1968 * Inspired by "How to Print Floating-Point Numbers Accurately" by
1969 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
1971 * Modifications:
1972 * 1. Rather than iterating, we use a simple numeric overestimate
1973 * to determine k = floor(log10(d)). We scale relevant
1974 * quantities using O(log2(k)) rather than O(k) multiplications.
1975 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
1976 * try to generate digits strictly left to right. Instead, we
1977 * compute with fewer bits and propagate the carry if necessary
1978 * when rounding the final digit up. This is often faster.
1979 * 3. Under the assumption that input will be rounded nearest,
1980 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
1981 * That is, we allow equality in stopping tests when the
1982 * round-nearest rule will give the same floating-point value
1983 * as would satisfaction of the stopping test with strict
1984 * inequality.
1985 * 4. We remove common factors of powers of 2 from relevant
1986 * quantities.
1987 * 5. When converting floating-point integers less than 1e16,
1988 * we use floating-point arithmetic rather than resorting
1989 * to multiple-precision integers.
1990 * 6. When asked to produce fewer than 15 digits, we first try
1991 * to get by with floating-point arithmetic; we resort to
1992 * multiple-precision integer arithmetic only if we cannot
1993 * guarantee that the floating-point calculation has given
1994 * the correctly rounded result. For k requested digits and
1995 * "uniformly" distributed input, the probability is
1996 * something like 10^(k-15) that we must resort to the Long
1997 * calculation.
2000 /* Additional notes (METD): (1) returns NULL on failure. (2) to avoid memory
2001 leakage, a successful call to _Py_dg_dtoa should always be matched by a
2002 call to _Py_dg_freedtoa. */
2004 char *
2005 _Py_dg_dtoa(double dd, int mode, int ndigits,
2006 int *decpt, int *sign, char **rve)
2008 /* Arguments ndigits, decpt, sign are similar to those
2009 of ecvt and fcvt; trailing zeros are suppressed from
2010 the returned string. If not null, *rve is set to point
2011 to the end of the return value. If d is +-Infinity or NaN,
2012 then *decpt is set to 9999.
2014 mode:
2015 0 ==> shortest string that yields d when read in
2016 and rounded to nearest.
2017 1 ==> like 0, but with Steele & White stopping rule;
2018 e.g. with IEEE P754 arithmetic , mode 0 gives
2019 1e23 whereas mode 1 gives 9.999999999999999e22.
2020 2 ==> max(1,ndigits) significant digits. This gives a
2021 return value similar to that of ecvt, except
2022 that trailing zeros are suppressed.
2023 3 ==> through ndigits past the decimal point. This
2024 gives a return value similar to that from fcvt,
2025 except that trailing zeros are suppressed, and
2026 ndigits can be negative.
2027 4,5 ==> similar to 2 and 3, respectively, but (in
2028 round-nearest mode) with the tests of mode 0 to
2029 possibly return a shorter string that rounds to d.
2030 With IEEE arithmetic and compilation with
2031 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
2032 as modes 2 and 3 when FLT_ROUNDS != 1.
2033 6-9 ==> Debugging modes similar to mode - 4: don't try
2034 fast floating-point estimate (if applicable).
2036 Values of mode other than 0-9 are treated as mode 0.
2038 Sufficient space is allocated to the return value
2039 to hold the suppressed trailing zeros.
2042 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
2043 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
2044 spec_case, try_quick;
2045 Long L;
2046 int denorm;
2047 ULong x;
2048 Bigint *b, *b1, *delta, *mlo, *mhi, *S;
2049 U d2, eps, u;
2050 double ds;
2051 char *s, *s0;
2053 /* set pointers to NULL, to silence gcc compiler warnings and make
2054 cleanup easier on error */
2055 mlo = mhi = b = S = 0;
2056 s0 = 0;
2058 u.d = dd;
2059 if (word0(&u) & Sign_bit) {
2060 /* set sign for everything, including 0's and NaNs */
2061 *sign = 1;
2062 word0(&u) &= ~Sign_bit; /* clear sign bit */
2064 else
2065 *sign = 0;
2067 /* quick return for Infinities, NaNs and zeros */
2068 if ((word0(&u) & Exp_mask) == Exp_mask)
2070 /* Infinity or NaN */
2071 *decpt = 9999;
2072 if (!word1(&u) && !(word0(&u) & 0xfffff))
2073 return nrv_alloc("Infinity", rve, 8);
2074 return nrv_alloc("NaN", rve, 3);
2076 if (!dval(&u)) {
2077 *decpt = 1;
2078 return nrv_alloc("0", rve, 1);
2081 /* compute k = floor(log10(d)). The computation may leave k
2082 one too large, but should never leave k too small. */
2083 b = d2b(&u, &be, &bbits);
2084 if (b == NULL)
2085 goto failed_malloc;
2086 if ((i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1)))) {
2087 dval(&d2) = dval(&u);
2088 word0(&d2) &= Frac_mask1;
2089 word0(&d2) |= Exp_11;
2091 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
2092 * log10(x) = log(x) / log(10)
2093 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
2094 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
2096 * This suggests computing an approximation k to log10(d) by
2098 * k = (i - Bias)*0.301029995663981
2099 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
2101 * We want k to be too large rather than too small.
2102 * The error in the first-order Taylor series approximation
2103 * is in our favor, so we just round up the constant enough
2104 * to compensate for any error in the multiplication of
2105 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
2106 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
2107 * adding 1e-13 to the constant term more than suffices.
2108 * Hence we adjust the constant term to 0.1760912590558.
2109 * (We could get a more accurate k by invoking log10,
2110 * but this is probably not worthwhile.)
2113 i -= Bias;
2114 denorm = 0;
2116 else {
2117 /* d is denormalized */
2119 i = bbits + be + (Bias + (P-1) - 1);
2120 x = i > 32 ? word0(&u) << (64 - i) | word1(&u) >> (i - 32)
2121 : word1(&u) << (32 - i);
2122 dval(&d2) = x;
2123 word0(&d2) -= 31*Exp_msk1; /* adjust exponent */
2124 i -= (Bias + (P-1) - 1) + 1;
2125 denorm = 1;
2127 ds = (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 +
2128 i*0.301029995663981;
2129 k = (int)ds;
2130 if (ds < 0. && ds != k)
2131 k--; /* want k = floor(ds) */
2132 k_check = 1;
2133 if (k >= 0 && k <= Ten_pmax) {
2134 if (dval(&u) < tens[k])
2135 k--;
2136 k_check = 0;
2138 j = bbits - i - 1;
2139 if (j >= 0) {
2140 b2 = 0;
2141 s2 = j;
2143 else {
2144 b2 = -j;
2145 s2 = 0;
2147 if (k >= 0) {
2148 b5 = 0;
2149 s5 = k;
2150 s2 += k;
2152 else {
2153 b2 -= k;
2154 b5 = -k;
2155 s5 = 0;
2157 if (mode < 0 || mode > 9)
2158 mode = 0;
2160 try_quick = 1;
2162 if (mode > 5) {
2163 mode -= 4;
2164 try_quick = 0;
2166 leftright = 1;
2167 ilim = ilim1 = -1; /* Values for cases 0 and 1; done here to */
2168 /* silence erroneous "gcc -Wall" warning. */
2169 switch(mode) {
2170 case 0:
2171 case 1:
2172 i = 18;
2173 ndigits = 0;
2174 break;
2175 case 2:
2176 leftright = 0;
2177 /* no break */
2178 case 4:
2179 if (ndigits <= 0)
2180 ndigits = 1;
2181 ilim = ilim1 = i = ndigits;
2182 break;
2183 case 3:
2184 leftright = 0;
2185 /* no break */
2186 case 5:
2187 i = ndigits + k + 1;
2188 ilim = i;
2189 ilim1 = i - 1;
2190 if (i <= 0)
2191 i = 1;
2193 s0 = rv_alloc(i);
2194 if (s0 == NULL)
2195 goto failed_malloc;
2196 s = s0;
2199 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
2201 /* Try to get by with floating-point arithmetic. */
2203 i = 0;
2204 dval(&d2) = dval(&u);
2205 k0 = k;
2206 ilim0 = ilim;
2207 ieps = 2; /* conservative */
2208 if (k > 0) {
2209 ds = tens[k&0xf];
2210 j = k >> 4;
2211 if (j & Bletch) {
2212 /* prevent overflows */
2213 j &= Bletch - 1;
2214 dval(&u) /= bigtens[n_bigtens-1];
2215 ieps++;
2217 for(; j; j >>= 1, i++)
2218 if (j & 1) {
2219 ieps++;
2220 ds *= bigtens[i];
2222 dval(&u) /= ds;
2224 else if ((j1 = -k)) {
2225 dval(&u) *= tens[j1 & 0xf];
2226 for(j = j1 >> 4; j; j >>= 1, i++)
2227 if (j & 1) {
2228 ieps++;
2229 dval(&u) *= bigtens[i];
2232 if (k_check && dval(&u) < 1. && ilim > 0) {
2233 if (ilim1 <= 0)
2234 goto fast_failed;
2235 ilim = ilim1;
2236 k--;
2237 dval(&u) *= 10.;
2238 ieps++;
2240 dval(&eps) = ieps*dval(&u) + 7.;
2241 word0(&eps) -= (P-1)*Exp_msk1;
2242 if (ilim == 0) {
2243 S = mhi = 0;
2244 dval(&u) -= 5.;
2245 if (dval(&u) > dval(&eps))
2246 goto one_digit;
2247 if (dval(&u) < -dval(&eps))
2248 goto no_digits;
2249 goto fast_failed;
2251 if (leftright) {
2252 /* Use Steele & White method of only
2253 * generating digits needed.
2255 dval(&eps) = 0.5/tens[ilim-1] - dval(&eps);
2256 for(i = 0;;) {
2257 L = (Long)dval(&u);
2258 dval(&u) -= L;
2259 *s++ = '0' + (int)L;
2260 if (dval(&u) < dval(&eps))
2261 goto ret1;
2262 if (1. - dval(&u) < dval(&eps))
2263 goto bump_up;
2264 if (++i >= ilim)
2265 break;
2266 dval(&eps) *= 10.;
2267 dval(&u) *= 10.;
2270 else {
2271 /* Generate ilim digits, then fix them up. */
2272 dval(&eps) *= tens[ilim-1];
2273 for(i = 1;; i++, dval(&u) *= 10.) {
2274 L = (Long)(dval(&u));
2275 if (!(dval(&u) -= L))
2276 ilim = i;
2277 *s++ = '0' + (int)L;
2278 if (i == ilim) {
2279 if (dval(&u) > 0.5 + dval(&eps))
2280 goto bump_up;
2281 else if (dval(&u) < 0.5 - dval(&eps)) {
2282 while(*--s == '0');
2283 s++;
2284 goto ret1;
2286 break;
2290 fast_failed:
2291 s = s0;
2292 dval(&u) = dval(&d2);
2293 k = k0;
2294 ilim = ilim0;
2297 /* Do we have a "small" integer? */
2299 if (be >= 0 && k <= Int_max) {
2300 /* Yes. */
2301 ds = tens[k];
2302 if (ndigits < 0 && ilim <= 0) {
2303 S = mhi = 0;
2304 if (ilim < 0 || dval(&u) <= 5*ds)
2305 goto no_digits;
2306 goto one_digit;
2308 for(i = 1;; i++, dval(&u) *= 10.) {
2309 L = (Long)(dval(&u) / ds);
2310 dval(&u) -= L*ds;
2311 *s++ = '0' + (int)L;
2312 if (!dval(&u)) {
2313 break;
2315 if (i == ilim) {
2316 dval(&u) += dval(&u);
2317 if (dval(&u) > ds || (dval(&u) == ds && L & 1)) {
2318 bump_up:
2319 while(*--s == '9')
2320 if (s == s0) {
2321 k++;
2322 *s = '0';
2323 break;
2325 ++*s++;
2327 break;
2330 goto ret1;
2333 m2 = b2;
2334 m5 = b5;
2335 if (leftright) {
2337 denorm ? be + (Bias + (P-1) - 1 + 1) :
2338 1 + P - bbits;
2339 b2 += i;
2340 s2 += i;
2341 mhi = i2b(1);
2342 if (mhi == NULL)
2343 goto failed_malloc;
2345 if (m2 > 0 && s2 > 0) {
2346 i = m2 < s2 ? m2 : s2;
2347 b2 -= i;
2348 m2 -= i;
2349 s2 -= i;
2351 if (b5 > 0) {
2352 if (leftright) {
2353 if (m5 > 0) {
2354 mhi = pow5mult(mhi, m5);
2355 if (mhi == NULL)
2356 goto failed_malloc;
2357 b1 = mult(mhi, b);
2358 Bfree(b);
2359 b = b1;
2360 if (b == NULL)
2361 goto failed_malloc;
2363 if ((j = b5 - m5)) {
2364 b = pow5mult(b, j);
2365 if (b == NULL)
2366 goto failed_malloc;
2369 else {
2370 b = pow5mult(b, b5);
2371 if (b == NULL)
2372 goto failed_malloc;
2375 S = i2b(1);
2376 if (S == NULL)
2377 goto failed_malloc;
2378 if (s5 > 0) {
2379 S = pow5mult(S, s5);
2380 if (S == NULL)
2381 goto failed_malloc;
2384 /* Check for special case that d is a normalized power of 2. */
2386 spec_case = 0;
2387 if ((mode < 2 || leftright)
2389 if (!word1(&u) && !(word0(&u) & Bndry_mask)
2390 && word0(&u) & (Exp_mask & ~Exp_msk1)
2392 /* The special case */
2393 b2 += Log2P;
2394 s2 += Log2P;
2395 spec_case = 1;
2399 /* Arrange for convenient computation of quotients:
2400 * shift left if necessary so divisor has 4 leading 0 bits.
2402 * Perhaps we should just compute leading 28 bits of S once
2403 * and for all and pass them and a shift to quorem, so it
2404 * can do shifts and ors to compute the numerator for q.
2406 if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f))
2407 i = 32 - i;
2408 #define iInc 28
2409 i = dshift(S, s2);
2410 b2 += i;
2411 m2 += i;
2412 s2 += i;
2413 if (b2 > 0) {
2414 b = lshift(b, b2);
2415 if (b == NULL)
2416 goto failed_malloc;
2418 if (s2 > 0) {
2419 S = lshift(S, s2);
2420 if (S == NULL)
2421 goto failed_malloc;
2423 if (k_check) {
2424 if (cmp(b,S) < 0) {
2425 k--;
2426 b = multadd(b, 10, 0); /* we botched the k estimate */
2427 if (b == NULL)
2428 goto failed_malloc;
2429 if (leftright) {
2430 mhi = multadd(mhi, 10, 0);
2431 if (mhi == NULL)
2432 goto failed_malloc;
2434 ilim = ilim1;
2437 if (ilim <= 0 && (mode == 3 || mode == 5)) {
2438 if (ilim < 0) {
2439 /* no digits, fcvt style */
2440 no_digits:
2441 k = -1 - ndigits;
2442 goto ret;
2444 else {
2445 S = multadd(S, 5, 0);
2446 if (S == NULL)
2447 goto failed_malloc;
2448 if (cmp(b, S) <= 0)
2449 goto no_digits;
2451 one_digit:
2452 *s++ = '1';
2453 k++;
2454 goto ret;
2456 if (leftright) {
2457 if (m2 > 0) {
2458 mhi = lshift(mhi, m2);
2459 if (mhi == NULL)
2460 goto failed_malloc;
2463 /* Compute mlo -- check for special case
2464 * that d is a normalized power of 2.
2467 mlo = mhi;
2468 if (spec_case) {
2469 mhi = Balloc(mhi->k);
2470 if (mhi == NULL)
2471 goto failed_malloc;
2472 Bcopy(mhi, mlo);
2473 mhi = lshift(mhi, Log2P);
2474 if (mhi == NULL)
2475 goto failed_malloc;
2478 for(i = 1;;i++) {
2479 dig = quorem(b,S) + '0';
2480 /* Do we yet have the shortest decimal string
2481 * that will round to d?
2483 j = cmp(b, mlo);
2484 delta = diff(S, mhi);
2485 if (delta == NULL)
2486 goto failed_malloc;
2487 j1 = delta->sign ? 1 : cmp(b, delta);
2488 Bfree(delta);
2489 if (j1 == 0 && mode != 1 && !(word1(&u) & 1)
2491 if (dig == '9')
2492 goto round_9_up;
2493 if (j > 0)
2494 dig++;
2495 *s++ = dig;
2496 goto ret;
2498 if (j < 0 || (j == 0 && mode != 1
2499 && !(word1(&u) & 1)
2500 )) {
2501 if (!b->x[0] && b->wds <= 1) {
2502 goto accept_dig;
2504 if (j1 > 0) {
2505 b = lshift(b, 1);
2506 if (b == NULL)
2507 goto failed_malloc;
2508 j1 = cmp(b, S);
2509 if ((j1 > 0 || (j1 == 0 && dig & 1))
2510 && dig++ == '9')
2511 goto round_9_up;
2513 accept_dig:
2514 *s++ = dig;
2515 goto ret;
2517 if (j1 > 0) {
2518 if (dig == '9') { /* possible if i == 1 */
2519 round_9_up:
2520 *s++ = '9';
2521 goto roundoff;
2523 *s++ = dig + 1;
2524 goto ret;
2526 *s++ = dig;
2527 if (i == ilim)
2528 break;
2529 b = multadd(b, 10, 0);
2530 if (b == NULL)
2531 goto failed_malloc;
2532 if (mlo == mhi) {
2533 mlo = mhi = multadd(mhi, 10, 0);
2534 if (mlo == NULL)
2535 goto failed_malloc;
2537 else {
2538 mlo = multadd(mlo, 10, 0);
2539 if (mlo == NULL)
2540 goto failed_malloc;
2541 mhi = multadd(mhi, 10, 0);
2542 if (mhi == NULL)
2543 goto failed_malloc;
2547 else
2548 for(i = 1;; i++) {
2549 *s++ = dig = quorem(b,S) + '0';
2550 if (!b->x[0] && b->wds <= 1) {
2551 goto ret;
2553 if (i >= ilim)
2554 break;
2555 b = multadd(b, 10, 0);
2556 if (b == NULL)
2557 goto failed_malloc;
2560 /* Round off last digit */
2562 b = lshift(b, 1);
2563 if (b == NULL)
2564 goto failed_malloc;
2565 j = cmp(b, S);
2566 if (j > 0 || (j == 0 && dig & 1)) {
2567 roundoff:
2568 while(*--s == '9')
2569 if (s == s0) {
2570 k++;
2571 *s++ = '1';
2572 goto ret;
2574 ++*s++;
2576 else {
2577 while(*--s == '0');
2578 s++;
2580 ret:
2581 Bfree(S);
2582 if (mhi) {
2583 if (mlo && mlo != mhi)
2584 Bfree(mlo);
2585 Bfree(mhi);
2587 ret1:
2588 Bfree(b);
2589 *s = 0;
2590 *decpt = k + 1;
2591 if (rve)
2592 *rve = s;
2593 return s0;
2594 failed_malloc:
2595 if (S)
2596 Bfree(S);
2597 if (mlo && mlo != mhi)
2598 Bfree(mlo);
2599 if (mhi)
2600 Bfree(mhi);
2601 if (b)
2602 Bfree(b);
2603 if (s0)
2604 _Py_dg_freedtoa(s0);
2605 return NULL;
2607 #ifdef __cplusplus
2609 #endif
2611 #endif /* PY_NO_SHORT_FLOAT_REPR */