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