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