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