Branch data 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 : 17987558 : Balloc(int k)
361 : : {
362 : : int x;
363 : : Bigint *rv;
364 : : unsigned int len;
365 : :
366 [ + + + + ]: 17987558 : if (k <= Kmax && (rv = freelist[k]))
367 : 17983081 : freelist[k] = rv->next;
368 : : else {
369 : 4477 : x = 1 << k;
370 : 4477 : len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
371 : 4477 : /sizeof(double);
372 [ + + + + ]: 4477 : if (k <= Kmax && pmem_next - private_mem + len <= (Py_ssize_t)PRIVATE_mem) {
373 : 4437 : rv = (Bigint*)pmem_next;
374 : 4437 : pmem_next += len;
375 : : }
376 : : else {
377 : 40 : rv = (Bigint*)MALLOC(len*sizeof(double));
378 [ - + ]: 40 : if (rv == NULL)
379 : 0 : return NULL;
380 : : }
381 : 4477 : rv->k = k;
382 : 4477 : rv->maxwds = x;
383 : : }
384 : 17987558 : rv->sign = rv->wds = 0;
385 : 17987558 : return rv;
386 : : }
387 : :
388 : : /* Free a Bigint allocated with Balloc */
389 : :
390 : : static void
391 : 24540805 : Bfree(Bigint *v)
392 : : {
393 [ + + ]: 24540805 : if (v) {
394 [ + + ]: 17987214 : if (v->k > Kmax)
395 : 4 : FREE((void*)v);
396 : : else {
397 : 17987210 : v->next = freelist[v->k];
398 : 17987210 : freelist[v->k] = v;
399 : : }
400 : : }
401 : 24540805 : }
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 : 5457414 : 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 : 5457414 : wds = b->wds;
461 : 5457414 : x = b->x;
462 : 5457414 : i = 0;
463 : 5457414 : carry = a;
464 : : do {
465 : 32608411 : y = *x * (ULLong)m + carry;
466 : 32608411 : carry = y >> 32;
467 : 32608411 : *x++ = (ULong)(y & FFFFFFFF);
468 : : }
469 [ + + ]: 32608411 : while(++i < wds);
470 [ + + ]: 5457414 : if (carry) {
471 [ + + ]: 158800 : if (wds >= b->maxwds) {
472 : 11160 : b1 = Balloc(b->k+1);
473 [ - + ]: 11160 : if (b1 == NULL){
474 : 0 : Bfree(b);
475 : 0 : return NULL;
476 : : }
477 : 11160 : Bcopy(b1, b);
478 : 11160 : Bfree(b);
479 : 11160 : b = b1;
480 : : }
481 : 158800 : b->x[wds++] = (ULong)carry;
482 : 158800 : b->wds = wds;
483 : : }
484 : 5457414 : 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 : 27047 : s2b(const char *s, int nd0, int nd, ULong y9)
495 : : {
496 : : Bigint *b;
497 : : int i, k;
498 : : Long x, y;
499 : :
500 : 27047 : x = (nd + 8) / 9;
501 [ + + ]: 52536 : for(k = 0, y = 1; x > y; y <<= 1, k++) ;
502 : 27047 : b = Balloc(k);
503 [ - + ]: 27047 : if (b == NULL)
504 : 0 : return NULL;
505 : 27047 : b->x[0] = y9;
506 : 27047 : b->wds = 1;
507 : :
508 [ + + ]: 27047 : if (nd <= 9)
509 : 4807 : return b;
510 : :
511 : 22240 : s += 9;
512 [ + + ]: 131710 : for (i = 9; i < nd0; i++) {
513 : 109470 : b = multadd(b, 10, *s++ - '0');
514 [ - + ]: 109470 : if (b == NULL)
515 : 0 : return NULL;
516 : : }
517 : 22240 : s++;
518 [ + + ]: 104349 : for(; i < nd; i++) {
519 : 82109 : b = multadd(b, 10, *s++ - '0');
520 [ - + ]: 82109 : if (b == NULL)
521 : 0 : return NULL;
522 : : }
523 : 22240 : return b;
524 : : }
525 : :
526 : : /* count leading 0 bits in the 32-bit integer x. */
527 : :
528 : : static int
529 : 1403580 : hi0bits(ULong x)
530 : : {
531 : 1403580 : int k = 0;
532 : :
533 [ + + ]: 1403580 : if (!(x & 0xffff0000)) {
534 : 1365696 : k = 16;
535 : 1365696 : x <<= 16;
536 : : }
537 [ + + ]: 1403580 : if (!(x & 0xff000000)) {
538 : 1353975 : k += 8;
539 : 1353975 : x <<= 8;
540 : : }
541 [ + + ]: 1403580 : if (!(x & 0xf0000000)) {
542 : 1347793 : k += 4;
543 : 1347793 : x <<= 4;
544 : : }
545 [ + + ]: 1403580 : if (!(x & 0xc0000000)) {
546 : 1337999 : k += 2;
547 : 1337999 : x <<= 2;
548 : : }
549 [ + + ]: 1403580 : if (!(x & 0x80000000)) {
550 : 1340059 : k++;
551 [ - + ]: 1340059 : if (!(x & 0x40000000))
552 : 0 : return 32;
553 : : }
554 : 1403580 : 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 : 4749705 : lo0bits(ULong *y)
562 : : {
563 : : int k;
564 : 4749705 : ULong x = *y;
565 : :
566 [ + + ]: 4749705 : if (x & 7) {
567 [ + + ]: 156462 : if (x & 1)
568 : 79651 : return 0;
569 [ + + ]: 76811 : if (x & 2) {
570 : 44495 : *y = x >> 1;
571 : 44495 : return 1;
572 : : }
573 : 32316 : *y = x >> 2;
574 : 32316 : return 2;
575 : : }
576 : 4593243 : k = 0;
577 [ + + ]: 4593243 : if (!(x & 0xffff)) {
578 : 4563385 : k = 16;
579 : 4563385 : x >>= 16;
580 : : }
581 [ + + ]: 4593243 : if (!(x & 0xff)) {
582 : 5689 : k += 8;
583 : 5689 : x >>= 8;
584 : : }
585 [ + + ]: 4593243 : if (!(x & 0xf)) {
586 : 4469290 : k += 4;
587 : 4469290 : x >>= 4;
588 : : }
589 [ + + ]: 4593243 : if (!(x & 0x3)) {
590 : 126152 : k += 2;
591 : 126152 : x >>= 2;
592 : : }
593 [ + + ]: 4593243 : if (!(x & 1)) {
594 : 125144 : k++;
595 : 125144 : x >>= 1;
596 [ - + ]: 125144 : if (!x)
597 : 0 : return 32;
598 : : }
599 : 4593243 : *y = x;
600 : 4593243 : return k;
601 : : }
602 : :
603 : : /* convert a small nonnegative integer to a Bigint */
604 : :
605 : : static Bigint *
606 : 1544509 : i2b(int i)
607 : : {
608 : : Bigint *b;
609 : :
610 : 1544509 : b = Balloc(1);
611 [ - + ]: 1544509 : if (b == NULL)
612 : 0 : return NULL;
613 : 1544509 : b->x[0] = i;
614 : 1544509 : b->wds = 1;
615 : 1544509 : 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 : 1665399 : 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 [ + + + - : 1665399 : if ((!a->x[0] && a->wds == 1) || (!b->x[0] && b->wds == 1)) {
+ + + + ]
631 : 857 : c = Balloc(0);
632 [ - + ]: 857 : if (c == NULL)
633 : 0 : return NULL;
634 : 857 : c->wds = 1;
635 : 857 : c->x[0] = 0;
636 : 857 : return c;
637 : : }
638 : :
639 [ + + ]: 1664542 : if (a->wds < b->wds) {
640 : 1434723 : c = a;
641 : 1434723 : a = b;
642 : 1434723 : b = c;
643 : : }
644 : 1664542 : k = a->k;
645 : 1664542 : wa = a->wds;
646 : 1664542 : wb = b->wds;
647 : 1664542 : wc = wa + wb;
648 [ + + ]: 1664542 : if (wc > a->maxwds)
649 : 1381403 : k++;
650 : 1664542 : c = Balloc(k);
651 [ - + ]: 1664542 : if (c == NULL)
652 : 0 : return NULL;
653 [ + + ]: 8002278 : for(x = c->x, xa = x + wc; x < xa; x++)
654 : 6337736 : *x = 0;
655 : 1664542 : xa = a->x;
656 : 1664542 : xae = xa + wa;
657 : 1664542 : xb = b->x;
658 : 1664542 : xbe = xb + wb;
659 : 1664542 : xc0 = c->x;
660 [ + + ]: 3665211 : for(; xb < xbe; xc0++) {
661 [ + + ]: 2000669 : if ((y = *xb++)) {
662 : 1999630 : x = xa;
663 : 1999630 : xc = xc0;
664 : 1999630 : carry = 0;
665 : : do {
666 : 7487304 : z = *x++ * (ULLong)y + *xc + carry;
667 : 7487304 : carry = z >> 32;
668 : 7487304 : *xc++ = (ULong)(z & FFFFFFFF);
669 : : }
670 [ + + ]: 7487304 : while(x < xae);
671 : 1999630 : *xc = (ULong)carry;
672 : : }
673 : : }
674 [ + - + + ]: 3208941 : for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
675 : 1664542 : c->wds = wc;
676 : 1664542 : 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 : 1395717 : 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 [ + + ]: 1395717 : if ((i = k & 3)) {
697 : 460980 : b = multadd(b, p05[i-1], 0);
698 [ - + ]: 460980 : if (b == NULL)
699 : 0 : return NULL;
700 : : }
701 : :
702 [ + + ]: 1395717 : if (!(k >>= 2))
703 : 29126 : return b;
704 : 1366591 : p5 = p5s;
705 [ + + ]: 1366591 : if (!p5) {
706 : : /* first time */
707 : 91 : p5 = i2b(625);
708 [ - + ]: 91 : if (p5 == NULL) {
709 : 0 : Bfree(b);
710 : 0 : return NULL;
711 : : }
712 : 91 : p5s = p5;
713 : 91 : p5->next = 0;
714 : : }
715 : : for(;;) {
716 [ + + ]: 4280563 : if (k & 1) {
717 : 1574437 : b1 = mult(b, p5);
718 : 1574437 : Bfree(b);
719 : 1574437 : b = b1;
720 [ - + ]: 1574437 : if (b == NULL)
721 : 0 : return NULL;
722 : : }
723 [ + + ]: 4280563 : if (!(k >>= 1))
724 : 1366591 : break;
725 : 2913972 : p51 = p5->next;
726 [ + + ]: 2913972 : if (!p51) {
727 : 253 : p51 = mult(p5,p5);
728 [ - + ]: 253 : if (p51 == NULL) {
729 : 0 : Bfree(b);
730 : 0 : return NULL;
731 : : }
732 : 253 : p51->next = 0;
733 : 253 : p5->next = p51;
734 : : }
735 : 2913972 : p5 = p51;
736 : : }
737 : 1366591 : 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 : 2999012 : lshift(Bigint *b, int k)
798 : : {
799 : : int i, k1, n, n1;
800 : : Bigint *b1;
801 : : ULong *x, *x1, *xe, z;
802 : :
803 [ + - + + : 2999012 : if (!k || (!b->x[0] && b->wds == 1))
+ + ]
804 : 1384 : return b;
805 : :
806 : 2997628 : n = k >> 5;
807 : 2997628 : k1 = b->k;
808 : 2997628 : n1 = n + b->wds + 1;
809 [ + + ]: 4744987 : for(i = b->maxwds; n1 > i; i <<= 1)
810 : 1747359 : k1++;
811 : 2997628 : b1 = Balloc(k1);
812 [ - + ]: 2997628 : if (b1 == NULL) {
813 : 0 : Bfree(b);
814 : 0 : return NULL;
815 : : }
816 : 2997628 : x1 = b1->x;
817 [ + + ]: 5633166 : for(i = 0; i < n; i++)
818 : 2635538 : *x1++ = 0;
819 : 2997628 : x = b->x;
820 : 2997628 : xe = x + b->wds;
821 [ + + ]: 2997628 : if (k &= 0x1f) {
822 : 2994115 : k1 = 32 - k;
823 : 2994115 : z = 0;
824 : : do {
825 : 5941211 : *x1++ = *x << k | z;
826 : 5941211 : z = *x++ >> k1;
827 : : }
828 [ + + ]: 5941211 : while(x < xe);
829 [ + + ]: 2994115 : if ((*x1 = z))
830 : 59636 : ++n1;
831 : : }
832 : : else do
833 : 15854 : *x1++ = *x++;
834 [ + + ]: 15854 : while(x < xe);
835 : 2997628 : b1->wds = n1 - 1;
836 : 2997628 : Bfree(b);
837 : 2997628 : 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 : 10695761 : cmp(Bigint *a, Bigint *b)
845 : : {
846 : : ULong *xa, *xa0, *xb, *xb0;
847 : : int i, j;
848 : :
849 : 10695761 : i = a->wds;
850 : 10695761 : j = b->wds;
851 : : #ifdef DEBUG
852 : : if (i > 1 && !a->x[i-1])
853 : : Bug("cmp called with a->x[a->wds-1] == 0");
854 : : if (j > 1 && !b->x[j-1])
855 : : Bug("cmp called with b->x[b->wds-1] == 0");
856 : : #endif
857 [ + + ]: 10695761 : if (i -= j)
858 : 2078951 : return i;
859 : 8616810 : xa0 = a->x;
860 : 8616810 : xa = xa0 + j;
861 : 8616810 : xb0 = b->x;
862 : 8616810 : xb = xb0 + j;
863 : : for(;;) {
864 [ + + ]: 8725219 : if (*--xa != *--xb)
865 [ + + ]: 8588445 : return *xa < *xb ? -1 : 1;
866 [ + + ]: 136774 : if (xa <= xa0)
867 : 28365 : break;
868 : : }
869 : 28365 : 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 : 2147048 : 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 : 2147048 : i = cmp(a,b);
885 [ + + ]: 2147048 : if (!i) {
886 : 1468 : c = Balloc(0);
887 [ - + ]: 1468 : if (c == NULL)
888 : 0 : return NULL;
889 : 1468 : c->wds = 1;
890 : 1468 : c->x[0] = 0;
891 : 1468 : return c;
892 : : }
893 [ + + ]: 2145580 : if (i < 0) {
894 : 56178 : c = a;
895 : 56178 : a = b;
896 : 56178 : b = c;
897 : 56178 : i = 1;
898 : : }
899 : : else
900 : 2089402 : i = 0;
901 : 2145580 : c = Balloc(a->k);
902 [ - + ]: 2145580 : if (c == NULL)
903 : 0 : return NULL;
904 : 2145580 : c->sign = i;
905 : 2145580 : wa = a->wds;
906 : 2145580 : xa = a->x;
907 : 2145580 : xae = xa + wa;
908 : 2145580 : wb = b->wds;
909 : 2145580 : xb = b->x;
910 : 2145580 : xbe = xb + wb;
911 : 2145580 : xc = c->x;
912 : 2145580 : borrow = 0;
913 : : do {
914 : 12963790 : y = (ULLong)*xa++ - *xb++ - borrow;
915 : 12963790 : borrow = y >> 32 & (ULong)1;
916 : 12963790 : *xc++ = (ULong)(y & FFFFFFFF);
917 : : }
918 [ + + ]: 12963790 : while(xb < xbe);
919 [ + + ]: 3197905 : while(xa < xae) {
920 : 1052325 : y = *xa++ - borrow;
921 : 1052325 : borrow = y >> 32 & (ULong)1;
922 : 1052325 : *xc++ = (ULong)(y & FFFFFFFF);
923 : : }
924 [ + + ]: 2190958 : while(!*--xc)
925 : 45378 : wa--;
926 : 2145580 : c->wds = wa;
927 : 2145580 : 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 : 17224 : ulp(U *x)
935 : : {
936 : : Long L;
937 : : U u;
938 : :
939 : 17224 : L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
940 : 17224 : word0(&u) = L;
941 : 17224 : word1(&u) = 0;
942 : 17224 : return dval(&u);
943 : : }
944 : :
945 : : /* Convert a Bigint to a double plus an exponent */
946 : :
947 : : static double
948 : 31972 : b2d(Bigint *a, int *e)
949 : : {
950 : : ULong *xa, *xa0, w, y, z;
951 : : int k;
952 : : U d;
953 : :
954 : 31972 : xa0 = a->x;
955 : 31972 : xa = xa0 + a->wds;
956 : 31972 : y = *--xa;
957 : : #ifdef DEBUG
958 : : if (!y) Bug("zero y in b2d");
959 : : #endif
960 : 31972 : k = hi0bits(y);
961 : 31972 : *e = 32 - k;
962 [ + + ]: 31972 : if (k < Ebits) {
963 : 5905 : word0(&d) = Exp_1 | y >> (Ebits - k);
964 [ + + ]: 5905 : w = xa > xa0 ? *--xa : 0;
965 : 5905 : word1(&d) = y << ((32-Ebits) + k) | w >> (Ebits - k);
966 : 5905 : goto ret_d;
967 : : }
968 [ + + ]: 26067 : z = xa > xa0 ? *--xa : 0;
969 [ + + ]: 26067 : if (k -= Ebits) {
970 : 25344 : word0(&d) = Exp_1 | y << k | z >> (32 - k);
971 [ + + ]: 25344 : y = xa > xa0 ? *--xa : 0;
972 : 25344 : word1(&d) = z << k | y >> (32 - k);
973 : : }
974 : : else {
975 : 723 : word0(&d) = Exp_1 | y;
976 : 723 : word1(&d) = z;
977 : : }
978 : 31972 : ret_d:
979 : 31972 : 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 : 36968 : sd2b(U *d, int scale, int *e)
1005 : : {
1006 : : Bigint *b;
1007 : :
1008 : 36968 : b = Balloc(1);
1009 [ - + ]: 36968 : if (b == NULL)
1010 : 0 : return NULL;
1011 : :
1012 : : /* First construct b and e assuming that scale == 0. */
1013 : 36968 : b->wds = 2;
1014 : 36968 : b->x[0] = word1(d);
1015 : 36968 : b->x[1] = word0(d) & Frac_mask;
1016 : 36968 : *e = Etiny - 1 + (int)((word0(d) & Exp_mask) >> Exp_shift);
1017 [ + + ]: 36968 : if (*e < Etiny)
1018 : 1384 : *e = Etiny;
1019 : : else
1020 : 35584 : b->x[1] |= Exp_msk1;
1021 : :
1022 : : /* Now adjust for scale, provided that b != 0. */
1023 [ + + + + : 36968 : if (scale && (b->x[0] || b->x[1])) {
+ + ]
1024 : 7632 : *e -= scale;
1025 [ + + ]: 7632 : if (*e < Etiny) {
1026 : 5310 : scale = Etiny - *e;
1027 : 5310 : *e = Etiny;
1028 : : /* We can't shift more than P-1 bits without shifting out a 1. */
1029 : : assert(0 < scale && scale <= P - 1);
1030 [ + + ]: 5310 : if (scale >= 32) {
1031 : : /* The bits shifted out should all be zero. */
1032 : : assert(b->x[0] == 0);
1033 : 3226 : b->x[0] = b->x[1];
1034 : 3226 : b->x[1] = 0;
1035 : 3226 : scale -= 32;
1036 : : }
1037 [ + + ]: 5310 : if (scale) {
1038 : : /* The bits shifted out should all be zero. */
1039 : : assert(b->x[0] << (32 - scale) == 0);
1040 : 5294 : b->x[0] = (b->x[0] >> scale) | (b->x[1] << (32 - scale));
1041 : 5294 : b->x[1] >>= scale;
1042 : : }
1043 : : }
1044 : : }
1045 : : /* Ensure b is normalized. */
1046 [ + + ]: 36968 : if (!b->x[1])
1047 : 4977 : b->wds = 1;
1048 : :
1049 : 36968 : 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 : 4749705 : 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 : 4749705 : b = Balloc(1);
1070 [ - + ]: 4749705 : if (b == NULL)
1071 : 0 : return NULL;
1072 : 4749705 : x = b->x;
1073 : :
1074 : 4749705 : z = word0(d) & Frac_mask;
1075 : 4749705 : word0(d) &= 0x7fffffff; /* clear sign bit, which we ignore */
1076 [ + + ]: 4749705 : if ((de = (int)(word0(d) >> Exp_shift)))
1077 : 4748700 : z |= Exp_msk1;
1078 [ + + ]: 4749705 : if ((y = word1(d))) {
1079 [ + + ]: 183786 : if ((k = lo0bits(&y))) {
1080 : 104181 : x[0] = y | z << (32 - k);
1081 : 104181 : z >>= k;
1082 : : }
1083 : : else
1084 : 79605 : x[0] = y;
1085 : 183786 : i =
1086 [ + + ]: 183786 : b->wds = (x[1] = z) ? 2 : 1;
1087 : : }
1088 : : else {
1089 : 4565919 : k = lo0bits(&z);
1090 : 4565919 : x[0] = z;
1091 : 4565919 : i =
1092 : 4565919 : b->wds = 1;
1093 : 4565919 : k += 32;
1094 : : }
1095 [ + + ]: 4749705 : if (de) {
1096 : 4748700 : *e = de - Bias - (P-1) + k;
1097 : 4748700 : *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 : 4749705 : 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 : 15986 : ratio(Bigint *a, Bigint *b)
1111 : : {
1112 : : U da, db;
1113 : : int k, ka, kb;
1114 : :
1115 : 15986 : dval(&da) = b2d(a, &ka);
1116 : 15986 : dval(&db) = b2d(b, &kb);
1117 : 15986 : k = ka - kb + 32*(a->wds - b->wds);
1118 [ + + ]: 15986 : if (k > 0)
1119 : 12994 : word0(&da) += k*Exp_msk1;
1120 : : else {
1121 : 2992 : k = -k;
1122 : 2992 : word0(&db) += k*Exp_msk1;
1123 : : }
1124 : 15986 : 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 : 1370603 : dshift(Bigint *b, int p2)
1152 : : {
1153 : 1370603 : int rv = hi0bits(b->x[b->wds-1]) - 4;
1154 [ + + ]: 1370603 : if (p2 > 0)
1155 : 1325062 : rv -= p2;
1156 : 1370603 : 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 : 2951145 : quorem(Bigint *b, Bigint *S)
1165 : : {
1166 : : int n;
1167 : : ULong *bx, *bxe, q, *sx, *sxe;
1168 : : ULLong borrow, carry, y, ys;
1169 : :
1170 : 2951145 : n = S->wds;
1171 : : #ifdef DEBUG
1172 : : /*debug*/ if (b->wds > n)
1173 : : /*debug*/ Bug("oversize b in quorem");
1174 : : #endif
1175 [ + + ]: 2951145 : if (b->wds < n)
1176 : 3082 : return 0;
1177 : 2948063 : sx = S->x;
1178 : 2948063 : sxe = sx + --n;
1179 : 2948063 : bx = b->x;
1180 : 2948063 : bxe = bx + n;
1181 : 2948063 : q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
1182 : : #ifdef DEBUG
1183 : : /*debug*/ if (q > 9)
1184 : : /*debug*/ Bug("oversized quotient in quorem");
1185 : : #endif
1186 [ + + ]: 2948063 : if (q) {
1187 : 2664244 : borrow = 0;
1188 : 2664244 : carry = 0;
1189 : : do {
1190 : 18908126 : ys = *sx++ * (ULLong)q + carry;
1191 : 18908126 : carry = ys >> 32;
1192 : 18908126 : y = *bx - (ys & FFFFFFFF) - borrow;
1193 : 18908126 : borrow = y >> 32 & (ULong)1;
1194 : 18908126 : *bx++ = (ULong)(y & FFFFFFFF);
1195 : : }
1196 [ + + ]: 18908126 : while(sx <= sxe);
1197 [ + + ]: 2664244 : if (!*bxe) {
1198 : 2 : bx = b->x;
1199 [ + - - + ]: 2 : while(--bxe > bx && !*bxe)
1200 : 0 : --n;
1201 : 2 : b->wds = n;
1202 : : }
1203 : : }
1204 [ + + ]: 2948063 : if (cmp(b, S) >= 0) {
1205 : 28080 : q++;
1206 : 28080 : borrow = 0;
1207 : 28080 : carry = 0;
1208 : 28080 : bx = b->x;
1209 : 28080 : sx = S->x;
1210 : : do {
1211 : 88930 : ys = *sx++ + carry;
1212 : 88930 : carry = ys >> 32;
1213 : 88930 : y = *bx - (ys & FFFFFFFF) - borrow;
1214 : 88930 : borrow = y >> 32 & (ULong)1;
1215 : 88930 : *bx++ = (ULong)(y & FFFFFFFF);
1216 : : }
1217 [ + + ]: 88930 : while(sx <= sxe);
1218 : 28080 : bx = b->x;
1219 : 28080 : bxe = bx + n;
1220 [ + + ]: 28080 : if (!*bxe) {
1221 [ + + + + ]: 50767 : while(--bxe > bx && !*bxe)
1222 : 22712 : --n;
1223 : 28055 : b->wds = n;
1224 : : }
1225 : : }
1226 : 2948063 : 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 : 1434 : sulp(U *x, BCinfo *bc)
1237 : : {
1238 : : U u;
1239 : :
1240 [ + + + + ]: 1434 : if (bc->scale && 2*P + 1 > (int)((word0(x) & Exp_mask) >> Exp_shift)) {
1241 : : /* rv/2^bc->scale is subnormal */
1242 : 196 : word0(&u) = (P+2)*Exp_msk1;
1243 : 196 : word1(&u) = 0;
1244 : 196 : return u.d;
1245 : : }
1246 : : else {
1247 : : assert(word0(x) || word1(x)); /* x != 0.0 */
1248 : 1238 : 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 : 3257 : 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 : 3257 : nd = bc->nd;
1305 : 3257 : nd0 = bc->nd0;
1306 : 3257 : p5 = nd + bc->e0;
1307 : 3257 : b = sd2b(rv, bc->scale, &p2);
1308 [ - + ]: 3257 : 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 : 3257 : 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 : 3257 : b = lshift(b, 1);
1318 [ - + ]: 3257 : if (b == NULL)
1319 : 0 : return -1;
1320 : 3257 : b->x[0] |= 1;
1321 : 3257 : p2--;
1322 : :
1323 : 3257 : p2 -= p5;
1324 : 3257 : d = i2b(1);
1325 [ - + ]: 3257 : 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 [ + + ]: 3257 : if (p5 > 0) {
1333 : 1078 : d = pow5mult(d, p5);
1334 [ - + ]: 1078 : if (d == NULL) {
1335 : 0 : Bfree(b);
1336 : 0 : return -1;
1337 : : }
1338 : : }
1339 [ + + ]: 2179 : else if (p5 < 0) {
1340 : 1863 : b = pow5mult(b, -p5);
1341 [ - + ]: 1863 : if (b == NULL) {
1342 : 0 : Bfree(d);
1343 : 0 : return -1;
1344 : : }
1345 : : }
1346 [ + + ]: 3257 : if (p2 > 0) {
1347 : 709 : b2 = p2;
1348 : 709 : d2 = 0;
1349 : : }
1350 : : else {
1351 : 2548 : b2 = 0;
1352 : 2548 : d2 = -p2;
1353 : : }
1354 : 3257 : i = dshift(d, d2);
1355 [ + + ]: 3257 : if ((b2 += i) > 0) {
1356 : 3240 : b = lshift(b, b2);
1357 [ - + ]: 3240 : if (b == NULL) {
1358 : 0 : Bfree(d);
1359 : 0 : return -1;
1360 : : }
1361 : : }
1362 [ + + ]: 3257 : if ((d2 += i) > 0) {
1363 : 3240 : d = lshift(d, d2);
1364 [ - + ]: 3240 : 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 [ + + ]: 3257 : if (cmp(b, d) >= 0)
1374 : : /* b/d >= 1 */
1375 : 79 : dd = -1;
1376 : : else {
1377 : 3178 : i = 0;
1378 : : for(;;) {
1379 : 333273 : b = multadd(b, 10, 0);
1380 [ - + ]: 333273 : if (b == NULL) {
1381 : 0 : Bfree(d);
1382 : 0 : return -1;
1383 : : }
1384 [ + + ]: 333273 : dd = s0[i < nd0 ? i : i+1] - '0' - quorem(b, d);
1385 : 333273 : i++;
1386 : :
1387 [ + + ]: 333273 : if (dd)
1388 : 2234 : break;
1389 [ + + + + ]: 331039 : if (!b->x[0] && b->wds == 1) {
1390 : : /* b/d == 0 */
1391 : 943 : dd = i < nd;
1392 : 943 : break;
1393 : : }
1394 [ + + ]: 330096 : if (!(i < nd)) {
1395 : : /* b/d != 0, but digits of s0 exhausted */
1396 : 1 : dd = -1;
1397 : 1 : break;
1398 : : }
1399 : : }
1400 : : }
1401 : 3257 : Bfree(b);
1402 : 3257 : Bfree(d);
1403 [ + + + + : 3257 : if (dd > 0 || (dd == 0 && odd))
+ + ]
1404 : 793 : dval(rv) += sulp(rv, bc);
1405 : 3257 : 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 : 3812 : _Py_dg_stdnan(int sign)
1417 : : {
1418 : : U rv;
1419 : 3812 : word0(&rv) = NAN_WORD0;
1420 : 3812 : word1(&rv) = NAN_WORD1;
1421 [ + + ]: 3812 : if (sign)
1422 : 19 : word0(&rv) |= Sign_bit;
1423 : 3812 : 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 : 5775 : _Py_dg_infinity(int sign)
1431 : : {
1432 : : U rv;
1433 : 5775 : word0(&rv) = POSINF_WORD0;
1434 : 5775 : word1(&rv) = POSINF_WORD1;
1435 [ + + ]: 5775 : return sign ? -dval(&rv) : dval(&rv);
1436 : : }
1437 : :
1438 : : double
1439 : 1337123 : _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 : 1337123 : Bigint *bb = NULL, *bd = NULL, *bd0 = NULL, *bs = NULL, *delta = NULL;
1450 : : size_t ndigits, fraclen;
1451 : : double result;
1452 : :
1453 : 1337123 : dval(&rv) = 0.;
1454 : :
1455 : : /* Start parsing. */
1456 : 1337123 : c = *(s = s00);
1457 : :
1458 : : /* Parse optional sign, if present. */
1459 : 1337123 : sign = 0;
1460 [ + + + ]: 1337123 : switch (c) {
1461 : 14183 : case '-':
1462 : 14183 : sign = 1;
1463 : : /* fall through */
1464 : 17714 : case '+':
1465 : 17714 : c = *++s;
1466 : : }
1467 : :
1468 : : /* Skip leading zeros: lz is true iff there were leading zeros. */
1469 : 1337123 : s1 = s;
1470 [ + + ]: 2591489 : while (c == '0')
1471 : 1254366 : c = *++s;
1472 : 1337123 : 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 : 1337123 : s0 = s1 = s;
1479 [ + + + + ]: 4129748 : while ('0' <= c && c <= '9')
1480 : 2792625 : c = *++s;
1481 : 1337123 : ndigits = s - s1;
1482 : 1337123 : fraclen = 0;
1483 : :
1484 : : /* Parse decimal point and following digits. */
1485 [ + + ]: 1337123 : if (c == '.') {
1486 : 84629 : c = *++s;
1487 [ + + ]: 84629 : if (!ndigits) {
1488 : 32683 : s1 = s;
1489 [ + + ]: 111109 : while (c == '0')
1490 : 78426 : c = *++s;
1491 [ + + + + ]: 32683 : lz = lz || s != s1;
1492 : 32683 : fraclen += (s - s1);
1493 : 32683 : s0 = s;
1494 : : }
1495 : 84629 : s1 = s;
1496 [ + + + + ]: 439294 : while ('0' <= c && c <= '9')
1497 : 354665 : c = *++s;
1498 : 84629 : ndigits += s - s1;
1499 : 84629 : 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 [ + + + + ]: 1337123 : if (!ndigits && !lz) {
1506 [ + - ]: 7963 : if (se)
1507 : 7963 : *se = (char *)s00;
1508 : 7963 : 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 [ + - - + ]: 1329160 : if (ndigits > MAX_DIGITS || fraclen > MAX_DIGITS) {
1514 [ # # ]: 0 : if (se)
1515 : 0 : *se = (char *)s00;
1516 : 0 : goto parse_error;
1517 : : }
1518 : 1329160 : nd = (int)ndigits;
1519 : 1329160 : nd0 = (int)ndigits - (int)fraclen;
1520 : :
1521 : : /* Parse exponent. */
1522 : 1329160 : e = 0;
1523 [ + + + + ]: 1329160 : if (c == 'e' || c == 'E') {
1524 : 1246133 : s00 = s;
1525 : 1246133 : c = *++s;
1526 : :
1527 : : /* Exponent sign. */
1528 : 1246133 : esign = 0;
1529 [ + + + ]: 1246133 : switch (c) {
1530 : 1231830 : case '-':
1531 : 1231830 : esign = 1;
1532 : : /* fall through */
1533 : 1236464 : case '+':
1534 : 1236464 : c = *++s;
1535 : : }
1536 : :
1537 : : /* Skip zeros. lz is true iff there are leading zeros. */
1538 : 1246133 : s1 = s;
1539 [ + + ]: 1251580 : while (c == '0')
1540 : 5447 : c = *++s;
1541 : 1246133 : lz = s != s1;
1542 : :
1543 : : /* Get absolute value of the exponent. */
1544 : 1246133 : s1 = s;
1545 : 1246133 : abs_exp = 0;
1546 [ + + + + ]: 2524229 : while ('0' <= c && c <= '9') {
1547 : 1278096 : abs_exp = 10*abs_exp + (c - '0');
1548 : 1278096 : 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 [ + - - + ]: 1246133 : if (s - s1 > 9 || abs_exp > MAX_ABS_EXP)
1555 : 0 : e = (int)MAX_ABS_EXP;
1556 : : else
1557 : 1246133 : e = (int)abs_exp;
1558 [ + + ]: 1246133 : if (esign)
1559 : 1231830 : e = -e;
1560 : :
1561 : : /* A valid exponent must have at least one digit. */
1562 [ + + + + ]: 1246133 : if (s == s1 && !lz)
1563 : 1 : s = s00;
1564 : : }
1565 : :
1566 : : /* Adjust exponent to take into account position of the point. */
1567 : 1329160 : e -= nd - nd0;
1568 [ + + ]: 1329160 : if (nd0 <= 0)
1569 : 1246991 : nd0 = nd;
1570 : :
1571 : : /* Finished parsing. Set se to indicate how far we parsed */
1572 [ + + ]: 1329160 : if (se)
1573 : 112131 : *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 [ + + ]: 1329160 : if (!nd)
1578 : 1227792 : goto ret;
1579 [ + - ]: 313714 : for (i = nd; i > 0; ) {
1580 : 313714 : --i;
1581 [ + + + + ]: 313714 : if (s0[i < nd0 ? i : i+1] != '0') {
1582 : 101368 : ++i;
1583 : 101368 : break;
1584 : : }
1585 : : }
1586 : 101368 : e += nd - i;
1587 : 101368 : nd = i;
1588 [ + + ]: 101368 : if (nd0 > nd)
1589 : 9213 : 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 : 101368 : bc.e0 = e1 = e;
1628 : 101368 : y = z = 0;
1629 [ + + ]: 659700 : for (i = 0; i < nd; i++) {
1630 [ + + ]: 577301 : if (i < 9)
1631 [ + + ]: 393752 : y = 10*y + s0[i < nd0 ? i : i+1] - '0';
1632 [ + + ]: 183549 : else if (i < DBL_DIG+1)
1633 [ + + ]: 164580 : z = 10*z + s0[i < nd0 ? i : i+1] - '0';
1634 : : else
1635 : 18969 : break;
1636 : : }
1637 : :
1638 : 101368 : k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
1639 : 101368 : dval(&rv) = y;
1640 [ + + ]: 101368 : if (k > 9) {
1641 : 26521 : dval(&rv) = tens[k - 9] * dval(&rv) + z;
1642 : : }
1643 [ + + ]: 101368 : if (nd <= DBL_DIG
1644 : : && Flt_Rounds == 1
1645 : : ) {
1646 [ + + ]: 80325 : if (!e)
1647 : 24890 : goto ret;
1648 [ + + ]: 55435 : if (e > 0) {
1649 [ + + ]: 13774 : if (e <= Ten_pmax) {
1650 : 9257 : dval(&rv) *= tens[e];
1651 : 9257 : goto ret;
1652 : : }
1653 : 4517 : i = DBL_DIG - nd;
1654 [ + + ]: 4517 : if (e <= Ten_pmax + i) {
1655 : : /* A fancier test would sometimes let us do
1656 : : * this for larger i values.
1657 : : */
1658 : 903 : e -= i;
1659 : 903 : dval(&rv) *= tens[i];
1660 : 903 : dval(&rv) *= tens[e];
1661 : 903 : goto ret;
1662 : : }
1663 : : }
1664 [ + + ]: 41661 : else if (e >= -Ten_pmax) {
1665 : 38292 : dval(&rv) /= tens[-e];
1666 : 38292 : goto ret;
1667 : : }
1668 : : }
1669 : 28026 : e1 += nd - k;
1670 : :
1671 : 28026 : bc.scale = 0;
1672 : :
1673 : : /* Get starting approximation = rv * 10**e1 */
1674 : :
1675 [ + + ]: 28026 : if (e1 > 0) {
1676 [ + + ]: 8352 : if ((i = e1 & 15))
1677 : 8070 : dval(&rv) *= tens[i];
1678 [ + + ]: 8352 : if (e1 &= ~15) {
1679 [ + + ]: 7136 : if (e1 > DBL_MAX_10_EXP)
1680 : 770 : goto ovfl;
1681 : 6366 : e1 >>= 4;
1682 [ + + ]: 22037 : for(j = 0; e1 > 1; j++, e1 >>= 1)
1683 [ + + ]: 15671 : if (e1 & 1)
1684 : 5526 : dval(&rv) *= bigtens[j];
1685 : : /* The last multiplication could overflow. */
1686 : 6366 : word0(&rv) -= P*Exp_msk1;
1687 : 6366 : dval(&rv) *= bigtens[j];
1688 [ + + ]: 6366 : if ((z = word0(&rv) & Exp_mask)
1689 : : > Exp_msk1*(DBL_MAX_EXP+Bias-P))
1690 : 72 : goto ovfl;
1691 [ + + ]: 6294 : if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
1692 : : /* set to largest number */
1693 : : /* (Can't trust DBL_MAX) */
1694 : 425 : word0(&rv) = Big0;
1695 : 425 : word1(&rv) = Big1;
1696 : : }
1697 : : else
1698 : 5869 : word0(&rv) += P*Exp_msk1;
1699 : : }
1700 : : }
1701 [ + + ]: 19674 : 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 : 19280 : e1 = -e1;
1713 [ + + ]: 19280 : if ((i = e1 & 15))
1714 : 16367 : dval(&rv) /= tens[i];
1715 [ + + ]: 19280 : if (e1 >>= 4) {
1716 [ + + ]: 12094 : if (e1 >= 1 << n_bigtens)
1717 : 137 : goto undfl;
1718 [ + + ]: 11957 : if (e1 & Scale_Bit)
1719 : 4604 : bc.scale = 2*P;
1720 [ + + ]: 49713 : for(j = 0; e1 > 0; j++, e1 >>= 1)
1721 [ + + ]: 37756 : if (e1 & 1)
1722 : 22340 : dval(&rv) *= tinytens[j];
1723 [ + + ]: 11957 : if (bc.scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask)
1724 [ + + ]: 4604 : >> Exp_shift)) > 0) {
1725 : : /* scaled rv is denormal; clear j low bits */
1726 [ + + ]: 3562 : if (j >= 32) {
1727 : 2323 : word1(&rv) = 0;
1728 [ + + ]: 2323 : if (j >= 53)
1729 : 1232 : word0(&rv) = (P+2)*Exp_msk1;
1730 : : else
1731 : 1091 : word0(&rv) &= 0xffffffff << (j-32);
1732 : : }
1733 : : else
1734 : 1239 : word1(&rv) &= 0xffffffff << j;
1735 : : }
1736 [ - + ]: 11957 : 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 : 27047 : bc.nd = nd;
1746 : 27047 : 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 [ + + ]: 27047 : 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 [ + - ]: 7563 : 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 : 7563 : --i;
1761 [ + + + + ]: 7563 : if (s0[i < nd0 ? i : i+1] != '0') {
1762 : 6625 : ++i;
1763 : 6625 : break;
1764 : : }
1765 : : }
1766 : 6625 : e += nd - i;
1767 : 6625 : nd = i;
1768 [ + + ]: 6625 : if (nd0 > nd)
1769 : 5456 : nd0 = nd;
1770 [ + + ]: 6625 : if (nd < 9) { /* must recompute y */
1771 : 12 : y = 0;
1772 [ + + ]: 25 : for(i = 0; i < nd0; ++i)
1773 : 13 : y = 10*y + s0[i] - '0';
1774 [ + + ]: 13 : for(; i < nd; ++i)
1775 : 1 : y = 10*y + s0[i+1] - '0';
1776 : : }
1777 : : }
1778 : 27047 : bd0 = s2b(s0, nd0, nd, y);
1779 [ - + ]: 27047 : 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 : 33711 : bd = Balloc(bd0->k);
1811 [ - + ]: 33711 : if (bd == NULL) {
1812 : 0 : goto failed_malloc;
1813 : : }
1814 : 33711 : Bcopy(bd, bd0);
1815 : 33711 : bb = sd2b(&rv, bc.scale, &bbe); /* srv = bb * 2^bbe */
1816 [ - + ]: 33711 : 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 : 33711 : odd = bb->x[0] & 1;
1822 : :
1823 : : /* tdv = bd * 10**e; srv = bb * 2**bbe */
1824 : 33711 : bs = i2b(1);
1825 [ - + ]: 33711 : if (bs == NULL) {
1826 : 0 : goto failed_malloc;
1827 : : }
1828 : :
1829 [ + + ]: 33711 : if (e >= 0) {
1830 : 9263 : bb2 = bb5 = 0;
1831 : 9263 : bd2 = bd5 = e;
1832 : : }
1833 : : else {
1834 : 24448 : bb2 = bb5 = -e;
1835 : 24448 : bd2 = bd5 = 0;
1836 : : }
1837 [ + + ]: 33711 : if (bbe >= 0)
1838 : 9361 : bb2 += bbe;
1839 : : else
1840 : 24350 : bd2 -= bbe;
1841 : 33711 : bs2 = bb2;
1842 : 33711 : bb2++;
1843 : 33711 : 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 : 33711 : i = bb2 < bd2 ? bb2 : bd2;
1864 [ + + ]: 33711 : if (i > bs2)
1865 : 23749 : i = bs2;
1866 [ + + ]: 33711 : if (i > 0) {
1867 : 33671 : bb2 -= i;
1868 : 33671 : bd2 -= i;
1869 : 33671 : bs2 -= i;
1870 : : }
1871 : :
1872 : : /* Scale bb, bd, bs by the appropriate powers of 2 and 5. */
1873 [ + + ]: 33711 : if (bb5 > 0) {
1874 : 24448 : bs = pow5mult(bs, bb5);
1875 [ - + ]: 24448 : if (bs == NULL) {
1876 : 0 : goto failed_malloc;
1877 : : }
1878 : 24448 : Bigint *bb1 = mult(bs, bb);
1879 : 24448 : Bfree(bb);
1880 : 24448 : bb = bb1;
1881 [ - + ]: 24448 : if (bb == NULL) {
1882 : 0 : goto failed_malloc;
1883 : : }
1884 : : }
1885 [ + - ]: 33711 : if (bb2 > 0) {
1886 : 33711 : bb = lshift(bb, bb2);
1887 [ - + ]: 33711 : if (bb == NULL) {
1888 : 0 : goto failed_malloc;
1889 : : }
1890 : : }
1891 [ + + ]: 33711 : if (bd5 > 0) {
1892 : 8595 : bd = pow5mult(bd, bd5);
1893 [ - + ]: 8595 : if (bd == NULL) {
1894 : 0 : goto failed_malloc;
1895 : : }
1896 : : }
1897 [ + + ]: 33711 : if (bd2 > 0) {
1898 : 23749 : bd = lshift(bd, bd2);
1899 [ - + ]: 23749 : if (bd == NULL) {
1900 : 0 : goto failed_malloc;
1901 : : }
1902 : : }
1903 [ + + ]: 33711 : if (bs2 > 0) {
1904 : 8921 : bs = lshift(bs, bs2);
1905 [ - + ]: 8921 : 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 : 33711 : delta = diff(bb, bd);
1915 [ - + ]: 33711 : if (delta == NULL) {
1916 : 0 : goto failed_malloc;
1917 : : }
1918 : 33711 : dsign = delta->sign;
1919 : 33711 : delta->sign = 0;
1920 : 33711 : i = cmp(delta, bs);
1921 [ + + + + ]: 33711 : if (bc.nd > nd && i <= 0) {
1922 [ + + ]: 6054 : if (dsign)
1923 : 3106 : 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 [ + + + + ]: 2948 : 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 : 567 : 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 [ + + ]: 567 : if (j - bc.scale >= 2) {
1943 : 151 : dval(&rv) -= 0.5 * sulp(&rv, &bc);
1944 : 151 : break; /* Use bigcomp. */
1945 : : }
1946 : : }
1947 : :
1948 : : {
1949 : 2797 : bc.nd = nd;
1950 : 2797 : i = -1; /* Discarded digits make delta smaller. */
1951 : : }
1952 : : }
1953 : :
1954 [ + + ]: 30454 : if (i < 0) {
1955 : : /* Error is less than half an ulp -- check for
1956 : : * special case of mantissa a power of two.
1957 : : */
1958 [ + + + + : 12132 : if (dsign || word1(&rv) || word0(&rv) & Bndry_mask
+ + ]
1959 [ + + ]: 705 : || (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1
1960 : : ) {
1961 : : break;
1962 : : }
1963 [ + + + + ]: 72 : if (!delta->x[0] && delta->wds <= 1) {
1964 : : /* exact result */
1965 : 22 : break;
1966 : : }
1967 : 50 : delta = lshift(delta,Log2P);
1968 [ - + ]: 50 : if (delta == NULL) {
1969 : 0 : goto failed_malloc;
1970 : : }
1971 [ + + ]: 50 : if (cmp(delta, bs) > 0)
1972 : 16 : goto drop_down;
1973 : 34 : break;
1974 : : }
1975 [ + + ]: 18322 : if (i == 0) {
1976 : : /* exactly half-way between */
1977 [ + + ]: 2336 : if (dsign) {
1978 [ - + ]: 1521 : 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 [ - + - - ]: 815 : else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) {
1994 : 0 : drop_down:
1995 : : /* boundary case -- decrement exponent */
1996 [ - + ]: 16 : 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 : 16 : L = (word0(&rv) & Exp_mask) - Exp_msk1;
2010 : 16 : word0(&rv) = L | Bndry_mask1;
2011 : 16 : word1(&rv) = 0xffffffff;
2012 : 16 : break;
2013 : : }
2014 [ + + ]: 2336 : if (!odd)
2015 : 1846 : break;
2016 [ + + ]: 490 : if (dsign)
2017 : 467 : dval(&rv) += sulp(&rv, &bc);
2018 : : else {
2019 : 23 : dval(&rv) -= sulp(&rv, &bc);
2020 [ - + ]: 23 : if (!dval(&rv)) {
2021 [ # # ]: 0 : if (bc.nd >nd)
2022 : 0 : break;
2023 : 0 : goto undfl;
2024 : : }
2025 : : }
2026 : : /* dsign = 1 - dsign; */
2027 : 490 : break;
2028 : : }
2029 [ + + ]: 15986 : if ((aadj = ratio(delta, bs)) <= 2.) {
2030 [ + + ]: 6490 : if (dsign)
2031 : 3921 : aadj = aadj1 = 1.;
2032 [ + + - + ]: 2569 : else if (word1(&rv) || word0(&rv) & Bndry_mask) {
2033 [ + + - + ]: 1712 : if (word1(&rv) == Tiny1 && !word0(&rv)) {
2034 [ # # ]: 0 : if (bc.nd >nd)
2035 : 0 : break;
2036 : 0 : goto undfl;
2037 : : }
2038 : 1712 : aadj = 1.;
2039 : 1712 : aadj1 = -1.;
2040 : : }
2041 : : else {
2042 : : /* special case -- power of FLT_RADIX to be */
2043 : : /* rounded down... */
2044 : :
2045 [ - + ]: 857 : if (aadj < 2./FLT_RADIX)
2046 : 0 : aadj = 1./FLT_RADIX;
2047 : : else
2048 : 857 : aadj *= 0.5;
2049 : 857 : aadj1 = -aadj;
2050 : : }
2051 : : }
2052 : : else {
2053 : 9496 : aadj *= 0.5;
2054 [ + + ]: 9496 : aadj1 = dsign ? aadj : -aadj;
2055 : : if (Flt_Rounds == 0)
2056 : : aadj1 += 0.5;
2057 : : }
2058 : 15986 : y = word0(&rv) & Exp_mask;
2059 : :
2060 : : /* Check for overflow */
2061 : :
2062 [ + + ]: 15986 : if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
2063 : 1479 : dval(&rv0) = dval(&rv);
2064 : 1479 : word0(&rv) -= P*Exp_msk1;
2065 : 1479 : adj.d = aadj1 * ulp(&rv);
2066 : 1479 : dval(&rv) += adj.d;
2067 [ + + ]: 1479 : if ((word0(&rv) & Exp_mask) >=
2068 : : Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
2069 [ + - + + ]: 751 : if (word0(&rv0) == Big0 && word1(&rv0) == Big1) {
2070 : 588 : goto ovfl;
2071 : : }
2072 : 163 : word0(&rv) = Big0;
2073 : 163 : word1(&rv) = Big1;
2074 : 163 : goto cont;
2075 : : }
2076 : : else
2077 : 728 : word0(&rv) += P*Exp_msk1;
2078 : : }
2079 : : else {
2080 [ + + + + ]: 14507 : if (bc.scale && y <= 2*P*Exp_msk1) {
2081 [ + - ]: 2331 : if (aadj <= 0x7fffffff) {
2082 [ + + ]: 2331 : if ((z = (ULong)aadj) <= 0)
2083 : 734 : z = 1;
2084 : 2331 : aadj = z;
2085 [ + + ]: 2331 : aadj1 = dsign ? aadj : -aadj;
2086 : : }
2087 : 2331 : dval(&aadj2) = aadj1;
2088 : 2331 : word0(&aadj2) += (2*P+1)*Exp_msk1 - y;
2089 : 2331 : aadj1 = dval(&aadj2);
2090 : : }
2091 : 14507 : adj.d = aadj1 * ulp(&rv);
2092 : 14507 : dval(&rv) += adj.d;
2093 : : }
2094 : 15235 : z = word0(&rv) & Exp_mask;
2095 [ + + ]: 15235 : if (bc.nd == nd) {
2096 [ + + ]: 10928 : if (!bc.scale)
2097 [ + + ]: 9800 : if (y == z) {
2098 : : /* Can we stop now? */
2099 : 9769 : L = (Long)aadj;
2100 : 9769 : aadj -= L;
2101 : : /* The tolerances below are conservative. */
2102 [ + + + + : 9769 : if (dsign || word1(&rv) || word0(&rv) & Bndry_mask) {
- + ]
2103 [ + + + + ]: 9748 : if (aadj < .4999999 || aadj > .5000001)
2104 : : break;
2105 : : }
2106 [ - + ]: 21 : else if (aadj < .4999999/FLT_RADIX)
2107 : 21 : break;
2108 : : }
2109 : : }
2110 : 5466 : cont:
2111 : 6664 : Bfree(bb); bb = NULL;
2112 : 6664 : Bfree(bd); bd = NULL;
2113 : 6664 : Bfree(bs); bs = NULL;
2114 : 6664 : Bfree(delta); delta = NULL;
2115 : : }
2116 [ + + ]: 26459 : if (bc.nd > nd) {
2117 : 3257 : error = bigcomp(&rv, s0, &bc);
2118 [ - + ]: 3257 : if (error)
2119 : 0 : goto failed_malloc;
2120 : : }
2121 : :
2122 [ + + ]: 26459 : if (bc.scale) {
2123 : 4604 : word0(&rv0) = Exp_1 - 2*P*Exp_msk1;
2124 : 4604 : word1(&rv0) = 0;
2125 : 4604 : dval(&rv) *= dval(&rv0);
2126 : : }
2127 : :
2128 : 21855 : ret:
2129 [ + + ]: 1327593 : result = sign ? -dval(&rv) : dval(&rv);
2130 : 1327593 : goto done;
2131 : :
2132 : 7963 : parse_error:
2133 : 7963 : result = 0.0;
2134 : 7963 : 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 : 1430 : ovfl:
2146 : 1430 : errno = ERANGE;
2147 : : /* Can't trust HUGE_VAL */
2148 : 1430 : word0(&rv) = Exp_mask;
2149 : 1430 : word1(&rv) = 0;
2150 [ + + ]: 1430 : result = sign ? -dval(&rv) : dval(&rv);
2151 : 1430 : goto done;
2152 : :
2153 : 1337123 : done:
2154 : 1337123 : Bfree(bb);
2155 : 1337123 : Bfree(bd);
2156 : 1337123 : Bfree(bs);
2157 : 1337123 : Bfree(bd0);
2158 : 1337123 : Bfree(delta);
2159 : 1337123 : return result;
2160 : :
2161 : : }
2162 : :
2163 : : static char *
2164 : 4771685 : rv_alloc(int i)
2165 : : {
2166 : : int j, k, *r;
2167 : :
2168 : 4771685 : j = sizeof(ULong);
2169 : 4771685 : for(k = 0;
2170 [ + + ]: 4822906 : sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= (unsigned)i;
2171 : 51221 : j <<= 1)
2172 : 51221 : k++;
2173 : 4771685 : r = (int*)Balloc(k);
2174 [ - + ]: 4771685 : if (r == NULL)
2175 : 0 : return NULL;
2176 : 4771685 : *r = k;
2177 : 4771685 : return (char *)(r+1);
2178 : : }
2179 : :
2180 : : static char *
2181 : 21980 : nrv_alloc(const char *s, char **rve, int n)
2182 : : {
2183 : : char *rv, *t;
2184 : :
2185 : 21980 : rv = rv_alloc(n);
2186 [ - + ]: 21980 : if (rv == NULL)
2187 : 0 : return NULL;
2188 : 21980 : t = rv;
2189 [ + + ]: 144421 : while((*t = *s++)) t++;
2190 [ + - ]: 21980 : if (rve)
2191 : 21980 : *rve = t;
2192 : 21980 : 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 : 4771685 : _Py_dg_freedtoa(char *s)
2203 : : {
2204 : 4771685 : Bigint *b = (Bigint *)((int *)s - 1);
2205 : 4771685 : b->maxwds = 1 << (b->k = *(int*)b);
2206 : 4771685 : Bfree(b);
2207 : 4771685 : }
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 : 4771685 : _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 : 4771685 : mlo = mhi = S = 0;
2299 : 4771685 : s0 = 0;
2300 : :
2301 : 4771685 : u.d = dd;
2302 [ + + ]: 4771685 : if (word0(&u) & Sign_bit) {
2303 : : /* set sign for everything, including 0's and NaNs */
2304 : 3414858 : *sign = 1;
2305 : 3414858 : word0(&u) &= ~Sign_bit; /* clear sign bit */
2306 : : }
2307 : : else
2308 : 1356827 : *sign = 0;
2309 : :
2310 : : /* quick return for Infinities, NaNs and zeros */
2311 [ + + ]: 4771685 : if ((word0(&u) & Exp_mask) == Exp_mask)
2312 : : {
2313 : : /* Infinity or NaN */
2314 : 15448 : *decpt = 9999;
2315 [ + - + + ]: 15448 : if (!word1(&u) && !(word0(&u) & 0xfffff))
2316 : 13913 : return nrv_alloc("Infinity", rve, 8);
2317 : 1535 : return nrv_alloc("NaN", rve, 3);
2318 : : }
2319 [ + + ]: 4756237 : if (!dval(&u)) {
2320 : 6532 : *decpt = 1;
2321 : 6532 : 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 : 4749705 : b = d2b(&u, &be, &bbits);
2327 [ - + ]: 4749705 : if (b == NULL)
2328 : 0 : goto failed_malloc;
2329 [ + + ]: 4749705 : if ((i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1)))) {
2330 : 4748700 : dval(&d2) = dval(&u);
2331 : 4748700 : word0(&d2) &= Frac_mask1;
2332 : 4748700 : 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 : 4748700 : i -= Bias;
2357 : 4748700 : 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 : 4749705 : ds = (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 +
2371 : 4749705 : i*0.301029995663981;
2372 : 4749705 : k = (int)ds;
2373 [ + + + - ]: 4749705 : if (ds < 0. && ds != k)
2374 : 1293513 : k--; /* want k = floor(ds) */
2375 : 4749705 : k_check = 1;
2376 [ + + + + ]: 4749705 : if (k >= 0 && k <= Ten_pmax) {
2377 [ + + ]: 3411462 : if (dval(&u) < tens[k])
2378 : 1312 : k--;
2379 : 3411462 : k_check = 0;
2380 : : }
2381 : 4749705 : j = bbits - i - 1;
2382 [ + + ]: 4749705 : if (j >= 0) {
2383 : 4683537 : b2 = 0;
2384 : 4683537 : s2 = j;
2385 : : }
2386 : : else {
2387 : 66168 : b2 = -j;
2388 : 66168 : s2 = 0;
2389 : : }
2390 [ + + ]: 4749705 : if (k >= 0) {
2391 : 3455053 : b5 = 0;
2392 : 3455053 : s5 = k;
2393 : 3455053 : s2 += k;
2394 : : }
2395 : : else {
2396 : 1294652 : b2 -= k;
2397 : 1294652 : b5 = -k;
2398 : 1294652 : s5 = 0;
2399 : : }
2400 [ + - - + ]: 4749705 : if (mode < 0 || mode > 9)
2401 : 0 : mode = 0;
2402 : :
2403 : 4749705 : try_quick = 1;
2404 : :
2405 [ - + ]: 4749705 : if (mode > 5) {
2406 : 0 : mode -= 4;
2407 : 0 : try_quick = 0;
2408 : : }
2409 : 4749705 : leftright = 1;
2410 : 4749705 : ilim = ilim1 = -1; /* Values for cases 0 and 1; done here to */
2411 : : /* silence erroneous "gcc -Wall" warning. */
2412 [ + + - + : 4749705 : switch(mode) {
- - ]
2413 : 149783 : case 0:
2414 : : case 1:
2415 : 149783 : i = 18;
2416 : 149783 : ndigits = 0;
2417 : 149783 : break;
2418 : 3350894 : case 2:
2419 : 3350894 : leftright = 0;
2420 : : /* fall through */
2421 : 3350894 : case 4:
2422 [ - + ]: 3350894 : if (ndigits <= 0)
2423 : 0 : ndigits = 1;
2424 : 3350894 : ilim = ilim1 = i = ndigits;
2425 : 3350894 : break;
2426 : 1249028 : case 3:
2427 : 1249028 : leftright = 0;
2428 : : /* fall through */
2429 : 1249028 : case 5:
2430 : 1249028 : i = ndigits + k + 1;
2431 : 1249028 : ilim = i;
2432 : 1249028 : ilim1 = i - 1;
2433 [ + + ]: 1249028 : if (i <= 0)
2434 : 1215506 : i = 1;
2435 : : }
2436 : 4749705 : s0 = rv_alloc(i);
2437 [ - + ]: 4749705 : if (s0 == NULL)
2438 : 0 : goto failed_malloc;
2439 : 4749705 : s = s0;
2440 : :
2441 : :
2442 [ + + + + : 4749705 : if (ilim >= 0 && ilim <= Quick_max && try_quick) {
+ - ]
2443 : :
2444 : : /* Try to get by with floating-point arithmetic. */
2445 : :
2446 : 3373809 : i = 0;
2447 : 3373809 : dval(&d2) = dval(&u);
2448 : 3373809 : k0 = k;
2449 : 3373809 : ilim0 = ilim;
2450 : 3373809 : ieps = 2; /* conservative */
2451 [ + + ]: 3373809 : if (k > 0) {
2452 : 11887 : ds = tens[k&0xf];
2453 : 11887 : j = k >> 4;
2454 [ + + ]: 11887 : if (j & Bletch) {
2455 : : /* prevent overflows */
2456 : 2 : j &= Bletch - 1;
2457 : 2 : dval(&u) /= bigtens[n_bigtens-1];
2458 : 2 : ieps++;
2459 : : }
2460 [ + + ]: 12327 : for(; j; j >>= 1, i++)
2461 [ + + ]: 440 : if (j & 1) {
2462 : 284 : ieps++;
2463 : 284 : ds *= bigtens[i];
2464 : : }
2465 : 11887 : dval(&u) /= ds;
2466 : : }
2467 [ + + ]: 3361922 : else if ((j1 = -k)) {
2468 : 15536 : dval(&u) *= tens[j1 & 0xf];
2469 [ + + ]: 15930 : 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 [ + + + + : 3373809 : if (k_check && dval(&u) < 1. && ilim > 0) {
+ + ]
2476 [ + + ]: 16 : if (ilim1 <= 0)
2477 : 9 : goto fast_failed;
2478 : 7 : ilim = ilim1;
2479 : 7 : k--;
2480 : 7 : dval(&u) *= 10.;
2481 : 7 : ieps++;
2482 : : }
2483 : 3373800 : dval(&eps) = ieps*dval(&u) + 7.;
2484 : 3373800 : word0(&eps) -= (P-1)*Exp_msk1;
2485 [ + + ]: 3373800 : if (ilim == 0) {
2486 : 3218 : S = mhi = 0;
2487 : 3218 : dval(&u) -= 5.;
2488 [ + + ]: 3218 : if (dval(&u) > dval(&eps))
2489 : 882 : goto one_digit;
2490 [ + + ]: 2336 : if (dval(&u) < -dval(&eps))
2491 : 2329 : goto no_digits;
2492 : 7 : goto fast_failed;
2493 : : }
2494 [ - + ]: 3370582 : 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 : 3370582 : dval(&eps) *= tens[ilim-1];
2516 : 3440229 : for(i = 1;; i++, dval(&u) *= 10.) {
2517 : 3440229 : L = (Long)(dval(&u));
2518 [ + + ]: 3440229 : if (!(dval(&u) -= L))
2519 : 3343775 : ilim = i;
2520 : 3440229 : *s++ = '0' + (int)L;
2521 [ + + ]: 3440229 : if (i == ilim) {
2522 [ + + ]: 3370582 : if (dval(&u) > 0.5 + dval(&eps))
2523 : 11964 : goto bump_up;
2524 [ + + ]: 3358618 : else if (dval(&u) < 0.5 - dval(&eps)) {
2525 [ + + ]: 3359906 : while(*--s == '0');
2526 : 3357424 : s++;
2527 : 3357424 : goto ret1;
2528 : : }
2529 : 1194 : break;
2530 : : }
2531 : : }
2532 : : }
2533 : 1210 : fast_failed:
2534 : 1210 : s = s0;
2535 : 1210 : dval(&u) = dval(&d2);
2536 : 1210 : k = k0;
2537 : 1210 : ilim = ilim0;
2538 : : }
2539 : :
2540 : : /* Do we have a "small" integer? */
2541 : :
2542 [ + + + + ]: 1377106 : if (be >= 0 && k <= Int_max) {
2543 : : /* Yes. */
2544 : 9760 : ds = tens[k];
2545 [ + + - + ]: 9760 : 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 : 30813 : for(i = 1;; i++, dval(&u) *= 10.) {
2552 : 30813 : L = (Long)(dval(&u) / ds);
2553 : 30813 : dval(&u) -= L*ds;
2554 : 30813 : *s++ = '0' + (int)L;
2555 [ + + ]: 30813 : if (!dval(&u)) {
2556 : 9742 : break;
2557 : : }
2558 [ + + ]: 21071 : 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 [ + + ]: 13387 : while(*--s == '9')
2563 [ + + ]: 1491 : if (s == s0) {
2564 : 76 : k++;
2565 : 76 : *s = '0';
2566 : 76 : break;
2567 : : }
2568 : 11972 : ++*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 : 11982 : break;
2579 : : }
2580 : : }
2581 : 21724 : goto ret1;
2582 : : }
2583 : :
2584 : 1367346 : m2 = b2;
2585 : 1367346 : m5 = b5;
2586 [ + + ]: 1367346 : if (leftright) {
2587 : 140104 : i =
2588 [ + + ]: 140104 : denorm ? be + (Bias + (P-1) - 1 + 1) :
2589 : 139101 : 1 + P - bbits;
2590 : 140104 : b2 += i;
2591 : 140104 : s2 += i;
2592 : 140104 : mhi = i2b(1);
2593 [ - + ]: 140104 : if (mhi == NULL)
2594 : 0 : goto failed_malloc;
2595 : : }
2596 [ + + + - ]: 1367346 : if (m2 > 0 && s2 > 0) {
2597 : 1337830 : i = m2 < s2 ? m2 : s2;
2598 : 1337830 : b2 -= i;
2599 : 1337830 : m2 -= i;
2600 : 1337830 : s2 -= i;
2601 : : }
2602 [ + + ]: 1367346 : if (b5 > 0) {
2603 [ + + ]: 1279446 : if (leftright) {
2604 [ + - ]: 66261 : if (m5 > 0) {
2605 : 66261 : mhi = pow5mult(mhi, m5);
2606 [ - + ]: 66261 : if (mhi == NULL)
2607 : 0 : goto failed_malloc;
2608 : 66261 : b1 = mult(mhi, b);
2609 : 66261 : Bfree(b);
2610 : 66261 : b = b1;
2611 [ - + ]: 66261 : if (b == NULL)
2612 : 0 : goto failed_malloc;
2613 : : }
2614 [ - + ]: 66261 : 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 : 1213185 : b = pow5mult(b, b5);
2622 [ - + ]: 1213185 : if (b == NULL)
2623 : 0 : goto failed_malloc;
2624 : : }
2625 : : }
2626 : 1367346 : S = i2b(1);
2627 [ - + ]: 1367346 : if (S == NULL)
2628 : 0 : goto failed_malloc;
2629 [ + + ]: 1367346 : if (s5 > 0) {
2630 : 80287 : S = pow5mult(S, s5);
2631 [ - + ]: 80287 : 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 : 1367346 : spec_case = 0;
2638 [ + + - + ]: 1367346 : if ((mode < 2 || leftright)
2639 : : ) {
2640 [ + + + + ]: 140104 : if (!word1(&u) && !(word0(&u) & Bndry_mask)
2641 [ + + ]: 2700 : && word0(&u) & (Exp_mask & ~Exp_msk1)
2642 : : ) {
2643 : : /* The special case */
2644 : 2698 : b2 += Log2P;
2645 : 2698 : s2 += Log2P;
2646 : 2698 : 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 : 1367346 : i = dshift(S, s2);
2659 : 1367346 : b2 += i;
2660 : 1367346 : m2 += i;
2661 : 1367346 : s2 += i;
2662 [ + + ]: 1367346 : if (b2 > 0) {
2663 : 1367308 : b = lshift(b, b2);
2664 [ - + ]: 1367308 : if (b == NULL)
2665 : 0 : goto failed_malloc;
2666 : : }
2667 [ + + ]: 1367346 : if (s2 > 0) {
2668 : 1366124 : S = lshift(S, s2);
2669 [ - + ]: 1366124 : if (S == NULL)
2670 : 0 : goto failed_malloc;
2671 : : }
2672 [ + + ]: 1367346 : if (k_check) {
2673 [ + + ]: 1323158 : if (cmp(b,S) < 0) {
2674 : 1609 : k--;
2675 : 1609 : b = multadd(b, 10, 0); /* we botched the k estimate */
2676 [ - + ]: 1609 : if (b == NULL)
2677 : 0 : goto failed_malloc;
2678 [ + + ]: 1609 : if (leftright) {
2679 : 1575 : mhi = multadd(mhi, 10, 0);
2680 [ - + ]: 1575 : if (mhi == NULL)
2681 : 0 : goto failed_malloc;
2682 : : }
2683 : 1609 : ilim = ilim1;
2684 : : }
2685 : : }
2686 [ + + + + : 1367346 : if (ilim <= 0 && (mode == 3 || mode == 5)) {
- + ]
2687 [ + + ]: 1212304 : if (ilim < 0) {
2688 : : /* no digits, fcvt style */
2689 : 1212288 : no_digits:
2690 : 1214618 : k = -1 - ndigits;
2691 : 1214618 : goto ret;
2692 : : }
2693 : : else {
2694 : 16 : S = multadd(S, 5, 0);
2695 [ - + ]: 16 : if (S == NULL)
2696 : 0 : goto failed_malloc;
2697 [ + + ]: 16 : if (cmp(b, S) <= 0)
2698 : 1 : goto no_digits;
2699 : : }
2700 : 15 : one_digit:
2701 : 897 : *s++ = '1';
2702 : 897 : k++;
2703 : 897 : goto ret;
2704 : : }
2705 [ + + ]: 155042 : if (leftright) {
2706 [ + + ]: 140104 : if (m2 > 0) {
2707 : 138695 : mhi = lshift(mhi, m2);
2708 [ - + ]: 138695 : 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 : 140104 : mlo = mhi;
2717 [ + + ]: 140104 : if (spec_case) {
2718 : 2698 : mhi = Balloc(mhi->k);
2719 [ - + ]: 2698 : if (mhi == NULL)
2720 : 0 : goto failed_malloc;
2721 : 2698 : Bcopy(mhi, mlo);
2722 : 2698 : mhi = lshift(mhi, Log2P);
2723 [ - + ]: 2698 : if (mhi == NULL)
2724 : 0 : goto failed_malloc;
2725 : : }
2726 : :
2727 : 2113337 : for(i = 1;;i++) {
2728 : 2113337 : dig = quorem(b,S) + '0';
2729 : : /* Do we yet have the shortest decimal string
2730 : : * that will round to d?
2731 : : */
2732 : 2113337 : j = cmp(b, mlo);
2733 : 2113337 : delta = diff(S, mhi);
2734 [ - + ]: 2113337 : if (delta == NULL)
2735 : 0 : goto failed_malloc;
2736 [ + + ]: 2113337 : j1 = delta->sign ? 1 : cmp(b, delta);
2737 : 2113337 : Bfree(delta);
2738 [ + + + - : 2113337 : if (j1 == 0 && mode != 1 && !(word1(&u) & 1)
+ + ]
2739 : : ) {
2740 [ + + ]: 1448 : if (dig == '9')
2741 : 13 : goto round_9_up;
2742 [ + + ]: 1435 : if (j > 0)
2743 : 249 : dig++;
2744 : 1435 : *s++ = dig;
2745 : 1435 : goto ret;
2746 : : }
2747 [ + + + + : 2111889 : if (j < 0 || (j == 0 && mode != 1
+ - ]
2748 [ + + ]: 749 : && !(word1(&u) & 1)
2749 : : )) {
2750 [ + + + + ]: 97471 : if (!b->x[0] && b->wds <= 1) {
2751 : 7769 : goto accept_dig;
2752 : : }
2753 [ + + ]: 89702 : if (j1 > 0) {
2754 : 44699 : b = lshift(b, 1);
2755 [ - + ]: 44699 : if (b == NULL)
2756 : 0 : goto failed_malloc;
2757 : 44699 : j1 = cmp(b, S);
2758 [ + + + + : 44699 : if ((j1 > 0 || (j1 == 0 && dig & 1))
+ + ]
2759 [ + + ]: 22541 : && dig++ == '9')
2760 : 98 : goto round_9_up;
2761 : : }
2762 : 89604 : accept_dig:
2763 : 97373 : *s++ = dig;
2764 : 97373 : goto ret;
2765 : : }
2766 [ + + ]: 2014418 : if (j1 > 0) {
2767 [ + + ]: 41185 : if (dig == '9') { /* possible if i == 1 */
2768 : 738 : round_9_up:
2769 : 849 : *s++ = '9';
2770 : 849 : goto roundoff;
2771 : : }
2772 : 40447 : *s++ = dig + 1;
2773 : 40447 : goto ret;
2774 : : }
2775 : 1973233 : *s++ = dig;
2776 [ - + ]: 1973233 : if (i == ilim)
2777 : 0 : break;
2778 : 1973233 : b = multadd(b, 10, 0);
2779 [ - + ]: 1973233 : if (b == NULL)
2780 : 0 : goto failed_malloc;
2781 [ + + ]: 1973233 : if (mlo == mhi) {
2782 : 1940914 : mlo = mhi = multadd(mhi, 10, 0);
2783 [ - + ]: 1940914 : if (mlo == NULL)
2784 : 0 : goto failed_malloc;
2785 : : }
2786 : : else {
2787 : 32319 : mlo = multadd(mlo, 10, 0);
2788 [ - + ]: 32319 : if (mlo == NULL)
2789 : 0 : goto failed_malloc;
2790 : 32319 : mhi = multadd(mhi, 10, 0);
2791 [ - + ]: 32319 : if (mhi == NULL)
2792 : 0 : goto failed_malloc;
2793 : : }
2794 : : }
2795 : : }
2796 : : else
2797 : 504535 : for(i = 1;; i++) {
2798 : 504535 : *s++ = dig = quorem(b,S) + '0';
2799 [ + + + + ]: 504535 : if (!b->x[0] && b->wds <= 1) {
2800 : 11618 : goto ret;
2801 : : }
2802 [ + + ]: 492917 : if (i >= ilim)
2803 : 3320 : break;
2804 : 489597 : b = multadd(b, 10, 0);
2805 [ - + ]: 489597 : if (b == NULL)
2806 : 0 : goto failed_malloc;
2807 : : }
2808 : :
2809 : : /* Round off last digit */
2810 : :
2811 : 3320 : b = lshift(b, 1);
2812 [ - + ]: 3320 : if (b == NULL)
2813 : 0 : goto failed_malloc;
2814 : 3320 : j = cmp(b, S);
2815 [ + + + + : 3320 : if (j > 0 || (j == 0 && dig & 1)) {
+ + ]
2816 : 2148 : roundoff:
2817 [ + + ]: 3140 : while(*--s == '9')
2818 [ + + ]: 1000 : if (s == s0) {
2819 : 857 : k++;
2820 : 857 : *s++ = '1';
2821 : 857 : goto ret;
2822 : : }
2823 : 2140 : ++*s++;
2824 : : }
2825 : : else {
2826 [ + + ]: 1361 : while(*--s == '0');
2827 : 1172 : s++;
2828 : : }
2829 : 1370557 : ret:
2830 : 1370557 : Bfree(S);
2831 [ + + ]: 1370557 : if (mhi) {
2832 [ + - + + ]: 140104 : if (mlo && mlo != mhi)
2833 : 2698 : Bfree(mlo);
2834 : 140104 : Bfree(mhi);
2835 : : }
2836 : 1230453 : ret1:
2837 : 4749705 : Bfree(b);
2838 : 4749705 : *s = 0;
2839 : 4749705 : *decpt = k + 1;
2840 [ + - ]: 4749705 : if (rve)
2841 : 4749705 : *rve = s;
2842 : 4749705 : 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
|
