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