/home/mdboom/Work/builds/cpython/Python/dtoa.c
Line | Count | Source (jump to first uncovered line) |
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_mask11.35k | Exp_msk11.35k *(DBL_MAX_EXP+Bias1.35k -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 | Balloc(int k) |
361 | { |
362 | int x; |
363 | Bigint *rv; |
364 | unsigned int len; |
365 | |
366 | if (k <= Kmax && (rv = freelist[k])11.2M ) Branch (366:9): [True: 11.2M, False: 4]
Branch (366:22): [True: 11.2M, False: 49]
|
367 | freelist[k] = rv->next; |
368 | else { |
369 | x = 1 << k; |
370 | len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1) |
371 | /sizeof(double); |
372 | if (k <= Kmax && pmem_next - private_mem + len <= (Py_ssize_t)49 PRIVATE_mem49 ) { Branch (372:13): [True: 49, False: 4]
Branch (372:26): [True: 30, False: 19]
|
373 | rv = (Bigint*)pmem_next; |
374 | pmem_next += len; |
375 | } |
376 | else { |
377 | rv = (Bigint*)MALLOC(len*sizeof(double)); |
378 | if (rv == NULL) Branch (378:17): [True: 0, False: 23]
|
379 | return NULL; |
380 | } |
381 | rv->k = k; |
382 | rv->maxwds = x; |
383 | } |
384 | rv->sign = rv->wds = 0; |
385 | return rv; |
386 | } |
387 | |
388 | /* Free a Bigint allocated with Balloc */ |
389 | |
390 | static void |
391 | Bfree(Bigint *v) |
392 | { |
393 | if (v) { Branch (393:9): [True: 11.2M, False: 6.28M]
|
394 | if (v->k > Kmax) Branch (394:13): [True: 4, False: 11.2M]
|
395 | FREE((void*)v); |
396 | else { |
397 | v->next = freelist[v->k]; |
398 | freelist[v->k] = v; |
399 | } |
400 | } |
401 | } |
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 | 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 | wds = b->wds; |
461 | x = b->x; |
462 | i = 0; |
463 | carry = a; |
464 | do { |
465 | y = *x * (ULLong)m + carry; |
466 | carry = y >> 32; |
467 | *x++ = (ULong)(y & FFFFFFFF); |
468 | } |
469 | while(++i < wds); Branch (469:11): [True: 27.1M, False: 5.40M]
|
470 | if (carry) { Branch (470:9): [True: 155k, False: 5.24M]
|
471 | if (wds >= b->maxwds) { Branch (471:13): [True: 10.6k, False: 144k]
|
472 | b1 = Balloc(b->k+1); |
473 | if (b1 == NULL){ Branch (473:17): [True: 0, False: 10.6k]
|
474 | Bfree(b); |
475 | return NULL; |
476 | } |
477 | Bcopy(b1, b); |
478 | Bfree(b); |
479 | b = b1; |
480 | } |
481 | b->x[wds++] = (ULong)carry; |
482 | b->wds = wds; |
483 | } |
484 | 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 | s2b(const char *s, int nd0, int nd, ULong y9) |
495 | { |
496 | Bigint *b; |
497 | int i, k; |
498 | Long x, y; |
499 | |
500 | x = (nd + 8) / 9; |
501 | for(k = 0, y = 1; x > y; y <<= 1, k++22.2k ) ;22.2k Branch (501:23): [True: 22.2k, False: 23.8k]
|
502 | b = Balloc(k); |
503 | if (b == NULL) Branch (503:9): [True: 0, False: 23.8k]
|
504 | return NULL; |
505 | b->x[0] = y9; |
506 | b->wds = 1; |
507 | |
508 | if (nd <= 9) Branch (508:9): [True: 3.69k, False: 20.1k]
|
509 | return b; |
510 | |
511 | s += 9; |
512 | for (i = 9; i < nd0; i++106k ) { Branch (512:17): [True: 106k, False: 20.1k]
|
513 | b = multadd(b, 10, *s++ - '0'); |
514 | if (b == NULL) Branch (514:13): [True: 0, False: 106k]
|
515 | return NULL; |
516 | } |
517 | s++; |
518 | for(; i < nd; i++65.4k ) { Branch (518:11): [True: 65.4k, False: 20.1k]
|
519 | b = multadd(b, 10, *s++ - '0'); |
520 | if (b == NULL) Branch (520:13): [True: 0, False: 65.4k]
|
521 | return NULL; |
522 | } |
523 | return b; |
524 | } |
525 | |
526 | /* count leading 0 bits in the 32-bit integer x. */ |
527 | |
528 | static int |
529 | hi0bits(ULong x) |
530 | { |
531 | int k = 0; |
532 | |
533 | if (!(x & 0xffff0000)) { Branch (533:9): [True: 1.35M, False: 37.6k]
|
534 | k = 16; |
535 | x <<= 16; |
536 | } |
537 | if (!(x & 0xff000000)) { Branch (537:9): [True: 1.34M, False: 49.2k]
|
538 | k += 8; |
539 | x <<= 8; |
540 | } |
541 | if (!(x & 0xf0000000)) { Branch (541:9): [True: 1.34M, False: 53.3k]
|
542 | k += 4; |
543 | x <<= 4; |
544 | } |
545 | if (!(x & 0xc0000000)) { Branch (545:9): [True: 1.33M, False: 63.4k]
|
546 | k += 2; |
547 | x <<= 2; |
548 | } |
549 | if (!(x & 0x80000000)) { Branch (549:9): [True: 1.33M, False: 62.3k]
|
550 | k++; |
551 | if (!(x & 0x40000000)) Branch (551:13): [True: 0, False: 1.33M]
|
552 | return 32; |
553 | } |
554 | 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 | lo0bits(ULong *y) |
562 | { |
563 | int k; |
564 | ULong x = *y; |
565 | |
566 | if (x & 7) { Branch (566:9): [True: 154k, False: 1.24M]
|
567 | if (x & 1) Branch (567:13): [True: 78.3k, False: 76.3k]
|
568 | return 0; |
569 | if (x & 2) { Branch (569:13): [True: 44.8k, False: 31.4k]
|
570 | *y = x >> 1; |
571 | return 1; |
572 | } |
573 | *y = x >> 2; |
574 | return 2; |
575 | } |
576 | k = 0; |
577 | if (!(x & 0xffff)) { Branch (577:9): [True: 1.21M, False: 26.7k]
|
578 | k = 16; |
579 | x >>= 16; |
580 | } |
581 | if (!(x & 0xff)) { Branch (581:9): [True: 3.93k, False: 1.24M]
|
582 | k += 8; |
583 | x >>= 8; |
584 | } |
585 | if (!(x & 0xf)) { Branch (585:9): [True: 1.12M, False: 119k]
|
586 | k += 4; |
587 | x >>= 4; |
588 | } |
589 | if (!(x & 0x3)) { Branch (589:9): [True: 121k, False: 1.12M]
|
590 | k += 2; |
591 | x >>= 2; |
592 | } |
593 | if (!(x & 1)) { Branch (593:9): [True: 121k, False: 1.12M]
|
594 | k++; |
595 | x >>= 1; |
596 | if (!x) Branch (596:13): [True: 0, False: 121k]
|
597 | return 32; |
598 | } |
599 | *y = x; |
600 | return k; |
601 | } |
602 | |
603 | /* convert a small nonnegative integer to a Bigint */ |
604 | |
605 | static Bigint * |
606 | i2b(int i) |
607 | { |
608 | Bigint *b; |
609 | |
610 | b = Balloc(1); |
611 | if (b == NULL) Branch (611:9): [True: 0, False: 1.53M]
|
612 | return NULL; |
613 | b->x[0] = i; |
614 | b->wds = 1; |
615 | 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 | 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 | if ((!a->x[0] && a->wds == 12.07k ) || (!b->x[0] && b->wds == 11.90k )) { Branch (630:10): [True: 2.07k, False: 1.65M]
Branch (630:22): [True: 0, False: 2.07k]
Branch (630:39): [True: 1.90k, False: 1.65M]
Branch (630:51): [True: 862, False: 1.04k]
|
631 | c = Balloc(0); |
632 | if (c == NULL) Branch (632:13): [True: 0, False: 862]
|
633 | return NULL; |
634 | c->wds = 1; |
635 | c->x[0] = 0; |
636 | return c; |
637 | } |
638 | |
639 | if (a->wds < b->wds) { Branch (639:9): [True: 1.42M, False: 222k]
|
640 | c = a; |
641 | a = b; |
642 | b = c; |
643 | } |
644 | k = a->k; |
645 | wa = a->wds; |
646 | wb = b->wds; |
647 | wc = wa + wb; |
648 | if (wc > a->maxwds) Branch (648:9): [True: 1.37M, False: 277k]
|
649 | k++; |
650 | c = Balloc(k); |
651 | if (c == NULL) Branch (651:9): [True: 0, False: 1.65M]
|
652 | return NULL; |
653 | for(x = c->x, xa = x + wc; 1.65M x < xa; x++6.29M ) Branch (653:32): [True: 6.29M, False: 1.65M]
|
654 | *x = 0; |
655 | xa = a->x; |
656 | xae = xa + wa; |
657 | xb = b->x; |
658 | xbe = xb + wb; |
659 | xc0 = c->x; |
660 | for(; xb < xbe; xc0++1.98M ) { Branch (660:11): [True: 1.98M, False: 1.65M]
|
661 | if ((y = *xb++)) { Branch (661:13): [True: 1.98M, False: 1.04k]
|
662 | x = xa; |
663 | xc = xc0; |
664 | carry = 0; |
665 | do { |
666 | z = *x++ * (ULLong)y + *xc + carry; |
667 | carry = z >> 32; |
668 | *xc++ = (ULong)(z & FFFFFFFF); |
669 | } |
670 | while(x < xae); Branch (670:19): [True: 5.43M, False: 1.98M]
|
671 | *xc = (ULong)carry; |
672 | } |
673 | } |
674 | for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc1.53M ) ;1.53M Branch (674:36): [True: 3.18M, False: 0]
Branch (674:46): [True: 1.53M, False: 1.65M]
|
675 | c->wds = wc; |
676 | 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 | 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 | if ((i = k & 3)) { Branch (696:9): [True: 457k, False: 932k]
|
697 | b = multadd(b, p05[i-1], 0); |
698 | if (b == NULL) Branch (698:13): [True: 0, False: 457k]
|
699 | return NULL; |
700 | } |
701 | |
702 | if (!(k >>= 2)) Branch (702:9): [True: 27.6k, False: 1.36M]
|
703 | return b; |
704 | p5 = p5s; |
705 | if (!p5) { Branch (705:9): [True: 1, False: 1.36M]
|
706 | /* first time */ |
707 | p5 = i2b(625); |
708 | if (p5 == NULL) { Branch (708:13): [True: 0, False: 1]
|
709 | Bfree(b); |
710 | return NULL; |
711 | } |
712 | p5s = p5; |
713 | p5->next = 0; |
714 | } |
715 | for(;;)1.36M { |
716 | if (k & 1) { Branch (716:13): [True: 1.56M, False: 2.70M]
|
717 | b1 = mult(b, p5); |
718 | Bfree(b); |
719 | b = b1; |
720 | if (b == NULL) Branch (720:17): [True: 0, False: 1.56M]
|
721 | return NULL; |
722 | } |
723 | if (!(k >>= 1)) Branch (723:13): [True: 1.36M, False: 2.90M]
|
724 | break; |
725 | p51 = p5->next; |
726 | if (!p51) { Branch (726:13): [True: 6, False: 2.90M]
|
727 | p51 = mult(p5,p5); |
728 | if (p51 == NULL) { Branch (728:17): [True: 0, False: 6]
|
729 | Bfree(b); |
730 | return NULL; |
731 | } |
732 | p51->next = 0; |
733 | p5->next = p51; |
734 | } |
735 | p5 = p51; |
736 | } |
737 | 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 | lshift(Bigint *b, int k) |
798 | { |
799 | int i, k1, n, n1; |
800 | Bigint *b1; |
801 | ULong *x, *x1, *xe, z; |
802 | |
803 | if (!k || (!b->x[0] && b->wds == 115.8k )) Branch (803:9): [True: 0, False: 2.98M]
Branch (803:16): [True: 15.8k, False: 2.96M]
Branch (803:28): [True: 1.39k, False: 14.4k]
|
804 | return b; |
805 | |
806 | n = k >> 5; |
807 | k1 = b->k; |
808 | n1 = n + b->wds + 1; |
809 | for(i = b->maxwds; n1 > i; i <<= 11.73M ) Branch (809:24): [True: 1.73M, False: 2.97M]
|
810 | k1++; |
811 | b1 = Balloc(k1); |
812 | if (b1 == NULL) { Branch (812:9): [True: 0, False: 2.97M]
|
813 | Bfree(b); |
814 | return NULL; |
815 | } |
816 | x1 = b1->x; |
817 | for(i = 0; i < n; i++2.60M ) Branch (817:16): [True: 2.60M, False: 2.97M]
|
818 | *x1++ = 0; |
819 | x = b->x; |
820 | xe = x + b->wds; |
821 | if (k &= 0x1f) { Branch (821:9): [True: 2.97M, False: 3.28k]
|
822 | k1 = 32 - k; |
823 | z = 0; |
824 | do { |
825 | *x1++ = *x << k | z; |
826 | z = *x++ >> k1; |
827 | } |
828 | while(x < xe); Branch (828:15): [True: 2.94M, False: 2.97M]
|
829 | if ((*x1 = z)) Branch (829:13): [True: 57.7k, False: 2.91M]
|
830 | ++n1; |
831 | } |
832 | else do |
833 | *x1++ = *x++; |
834 | while(x < xe); Branch (834:15): [True: 12.6k, False: 3.28k]
|
835 | b1->wds = n1 - 1; |
836 | Bfree(b); |
837 | 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 | cmp(Bigint *a, Bigint *b) |
845 | { |
846 | ULong *xa, *xa0, *xb, *xb0; |
847 | int i, j; |
848 | |
849 | i = a->wds; |
850 | 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 | if (i -= j) Branch (857:9): [True: 2.05M, False: 8.56M]
|
858 | return i; |
859 | xa0 = a->x; |
860 | xa = xa0 + j; |
861 | xb0 = b->x; |
862 | xb = xb0 + j; |
863 | for(;;) { |
864 | if (*--xa != *--xb) Branch (864:13): [True: 8.53M, False: 127k]
|
865 | return *xa < *xb ? -15.07M : 13.45M ; Branch (865:20): [True: 5.07M, False: 3.45M]
|
866 | if (xa <= xa0) Branch (866:13): [True: 27.3k, False: 100k]
|
867 | break; |
868 | } |
869 | 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 | 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 | i = cmp(a,b); |
885 | if (!i) { Branch (885:9): [True: 1.41k, False: 2.12M]
|
886 | c = Balloc(0); |
887 | if (c == NULL) Branch (887:13): [True: 0, False: 1.41k]
|
888 | return NULL; |
889 | c->wds = 1; |
890 | c->x[0] = 0; |
891 | return c; |
892 | } |
893 | if (i < 0) { Branch (893:9): [True: 55.4k, False: 2.07M]
|
894 | c = a; |
895 | a = b; |
896 | b = c; |
897 | i = 1; |
898 | } |
899 | else |
900 | i = 0; |
901 | c = Balloc(a->k); |
902 | if (c == NULL) Branch (902:9): [True: 0, False: 2.12M]
|
903 | return NULL; |
904 | c->sign = i; |
905 | wa = a->wds; |
906 | xa = a->x; |
907 | xae = xa + wa; |
908 | wb = b->wds; |
909 | xb = b->x; |
910 | xbe = xb + wb; |
911 | xc = c->x; |
912 | borrow = 0; |
913 | do { |
914 | y = (ULLong)*xa++ - *xb++ - borrow; |
915 | borrow = y >> 32 & (ULong)1; |
916 | *xc++ = (ULong)(y & FFFFFFFF); |
917 | } |
918 | while(xb < xbe); Branch (918:11): [True: 10.9M, False: 2.12M]
|
919 | while(xa < xae) { Branch (919:11): [True: 1.04M, False: 2.12M]
|
920 | y = *xa++ - borrow; |
921 | borrow = y >> 32 & (ULong)1; |
922 | *xc++ = (ULong)(y & FFFFFFFF); |
923 | } |
924 | while(!*--xc) Branch (924:11): [True: 40.5k, False: 2.12M]
|
925 | wa--; |
926 | c->wds = wa; |
927 | 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 | ulp(U *x) |
935 | { |
936 | Long L; |
937 | U u; |
938 | |
939 | L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1; |
940 | word0(&u) = L; |
941 | word1(&u) = 0; |
942 | return dval(&u); |
943 | } |
944 | |
945 | /* Convert a Bigint to a double plus an exponent */ |
946 | |
947 | static double |
948 | b2d(Bigint *a, int *e) |
949 | { |
950 | ULong *xa, *xa0, w, y, z; |
951 | int k; |
952 | U d; |
953 | |
954 | xa0 = a->x; |
955 | xa = xa0 + a->wds; |
956 | y = *--xa; |
957 | #ifdef DEBUG |
958 | if (!y) Bug("zero y in b2d"); |
959 | #endif |
960 | k = hi0bits(y); |
961 | *e = 32 - k; |
962 | if (k < Ebits) { Branch (962:9): [True: 5.58k, False: 23.7k]
|
963 | word0(&d) = Exp_1 | y >> (Ebits - k); |
964 | w = xa > xa0 ? *--xa5.22k : 0366 ; Branch (964:13): [True: 5.22k, False: 366]
|
965 | word1(&d) = y << ((32-Ebits) + k) | w >> (Ebits - k); |
966 | goto ret_d; |
967 | } |
968 | z = xa > xa0 ? *--xa21.8k : 01.90k ; Branch (968:9): [True: 21.8k, False: 1.90k]
|
969 | if (k -= Ebits) { Branch (969:9): [True: 23.0k, False: 674]
|
970 | word0(&d) = Exp_1 | y << k | z >> (32 - k); |
971 | y = xa > xa0 ? *--xa9.58k : 013.4k ; Branch (971:13): [True: 9.58k, False: 13.4k]
|
972 | word1(&d) = z << k | y >> (32 - k); |
973 | } |
974 | else { |
975 | word0(&d) = Exp_1 | y; |
976 | word1(&d) = z; |
977 | } |
978 | ret_d: |
979 | 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 | sd2b(U *d, int scale, int *e) |
1005 | { |
1006 | Bigint *b; |
1007 | |
1008 | b = Balloc(1); |
1009 | if (b == NULL) Branch (1009:9): [True: 0, False: 33.7k]
|
1010 | return NULL; |
1011 | |
1012 | /* First construct b and e assuming that scale == 0. */ |
1013 | b->wds = 2; |
1014 | b->x[0] = word1(d); |
1015 | b->x[1] = word0(d) & Frac_mask; |
1016 | *e = Etiny - 1 + (int)((word0(d) & Exp_mask) >> Exp_shift); |
1017 | if (*e < Etiny) Branch (1017:9): [True: 1.39k, False: 32.3k]
|
1018 | *e = Etiny; |
1019 | else |
1020 | b->x[1] |= Exp_msk1; |
1021 | |
1022 | /* Now adjust for scale, provided that b != 0. */ |
1023 | if (scale && (8.86k b->x[0]8.86k || b->x[1]4.90k )) { Branch (1023:9): [True: 8.86k, False: 24.8k]
Branch (1023:19): [True: 3.96k, False: 4.90k]
Branch (1023:30): [True: 3.50k, False: 1.39k]
|
1024 | *e -= scale; |
1025 | if (*e < Etiny) { Branch (1025:13): [True: 5.26k, False: 2.20k]
|
1026 | scale = Etiny - *e; |
1027 | *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 | if (scale >= 32) { Branch (1030:17): [True: 3.16k, False: 2.10k]
|
1031 | /* The bits shifted out should all be zero. */ |
1032 | assert(b->x[0] == 0); |
1033 | b->x[0] = b->x[1]; |
1034 | b->x[1] = 0; |
1035 | scale -= 32; |
1036 | } |
1037 | if (scale) { Branch (1037:17): [True: 5.24k, False: 16]
|
1038 | /* The bits shifted out should all be zero. */ |
1039 | assert(b->x[0] << (32 - scale) == 0); |
1040 | b->x[0] = (b->x[0] >> scale) | (b->x[1] << (32 - scale)); |
1041 | b->x[1] >>= scale; |
1042 | } |
1043 | } |
1044 | } |
1045 | /* Ensure b is normalized. */ |
1046 | if (!b->x[1]) Branch (1046:9): [True: 4.89k, False: 28.8k]
|
1047 | b->wds = 1; |
1048 | |
1049 | 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 | 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 | b = Balloc(1); |
1070 | if (b == NULL) Branch (1070:9): [True: 0, False: 1.39M]
|
1071 | return NULL; |
1072 | x = b->x; |
1073 | |
1074 | z = word0(d) & Frac_mask; |
1075 | word0(d) &= 0x7fffffff; /* clear sign bit, which we ignore */ |
1076 | if ((de = (int)(word0(d) >> Exp_shift))) Branch (1076:9): [True: 1.39M, False: 991]
|
1077 | z |= Exp_msk1; |
1078 | if ((y = word1(d))) { Branch (1078:9): [True: 179k, False: 1.22M]
|
1079 | if ((k = lo0bits(&y))) { Branch (1079:13): [True: 101k, False: 78.2k]
|
1080 | x[0] = y | z << (32 - k); |
1081 | z >>= k; |
1082 | } |
1083 | else |
1084 | x[0] = y; |
1085 | i = |
1086 | b->wds = (x[1] = z) ? 2178k : 11.11k ; Branch (1086:22): [True: 178k, False: 1.11k]
|
1087 | } |
1088 | else { |
1089 | k = lo0bits(&z); |
1090 | x[0] = z; |
1091 | i = |
1092 | b->wds = 1; |
1093 | k += 32; |
1094 | } |
1095 | if (de) { Branch (1095:9): [True: 1.39M, False: 991]
|
1096 | *e = de - Bias - (P-1) + k; |
1097 | *bits = P - k; |
1098 | } |
1099 | else { |
1100 | *e = de - Bias - (P-1) + 1 + k; |
1101 | *bits = 32*i - hi0bits(x[i-1]); |
1102 | } |
1103 | 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 | ratio(Bigint *a, Bigint *b) |
1111 | { |
1112 | U da, db; |
1113 | int k, ka, kb; |
1114 | |
1115 | dval(&da) = b2d(a, &ka); |
1116 | dval(&db) = b2d(b, &kb); |
1117 | k = ka - kb + 32*(a->wds - b->wds); |
1118 | if (k > 0) Branch (1118:9): [True: 11.9k, False: 2.69k]
|
1119 | word0(&da) += k*Exp_msk1; |
1120 | else { |
1121 | k = -k; |
1122 | word0(&db) += k*Exp_msk1; |
1123 | } |
1124 | 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 | dshift(Bigint *b, int p2) |
1152 | { |
1153 | int rv = hi0bits(b->x[b->wds-1]) - 4; |
1154 | if (p2 > 0) Branch (1154:9): [True: 1.32M, False: 45.4k]
|
1155 | rv -= p2; |
1156 | 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 | quorem(Bigint *b, Bigint *S) |
1165 | { |
1166 | int n; |
1167 | ULong *bx, *bxe, q, *sx, *sxe; |
1168 | ULLong borrow, carry, y, ys; |
1169 | |
1170 | n = S->wds; |
1171 | #ifdef DEBUG |
1172 | /*debug*/ if (b->wds > n) |
1173 | /*debug*/ Bug("oversize b in quorem"); |
1174 | #endif |
1175 | if (b->wds < n) Branch (1175:9): [True: 2.99k, False: 2.92M]
|
1176 | return 0; |
1177 | sx = S->x; |
1178 | sxe = sx + --n; |
1179 | bx = b->x; |
1180 | bxe = bx + n; |
1181 | 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 | if (q) { Branch (1186:9): [True: 2.64M, False: 278k]
|
1187 | borrow = 0; |
1188 | carry = 0; |
1189 | do { |
1190 | ys = *sx++ * (ULLong)q + carry; |
1191 | carry = ys >> 32; |
1192 | y = *bx - (ys & FFFFFFFF) - borrow; |
1193 | borrow = y >> 32 & (ULong)1; |
1194 | *bx++ = (ULong)(y & FFFFFFFF); |
1195 | } |
1196 | while(sx <= sxe); Branch (1196:15): [True: 16.1M, False: 2.64M]
|
1197 | if (!*bxe) { Branch (1197:13): [True: 2, False: 2.64M]
|
1198 | bx = b->x; |
1199 | while(--bxe > bx && !*bxe) Branch (1199:19): [True: 2, False: 0]
Branch (1199:33): [True: 0, False: 2]
|
1200 | --n; |
1201 | b->wds = n; |
1202 | } |
1203 | } |
1204 | if (cmp(b, S) >= 0) { Branch (1204:9): [True: 25.9k, False: 2.89M]
|
1205 | q++; |
1206 | borrow = 0; |
1207 | carry = 0; |
1208 | bx = b->x; |
1209 | sx = S->x; |
1210 | do { |
1211 | ys = *sx++ + carry; |
1212 | carry = ys >> 32; |
1213 | y = *bx - (ys & FFFFFFFF) - borrow; |
1214 | borrow = y >> 32 & (ULong)1; |
1215 | *bx++ = (ULong)(y & FFFFFFFF); |
1216 | } |
1217 | while(sx <= sxe); Branch (1217:15): [True: 57.1k, False: 25.9k]
|
1218 | bx = b->x; |
1219 | bxe = bx + n; |
1220 | if (!*bxe) { Branch (1220:13): [True: 25.9k, False: 17]
|
1221 | while(--bxe > bx && !*bxe23.6k ) Branch (1221:19): [True: 23.6k, False: 24.8k]
Branch (1221:33): [True: 22.4k, False: 1.15k]
|
1222 | --n; |
1223 | b->wds = n; |
1224 | } |
1225 | } |
1226 | 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 | sulp(U *x, BCinfo *bc) |
1237 | { |
1238 | U u; |
1239 | |
1240 | if (bc->scale && 2*258 P258 + 1 > (int)((word0258 (x) & Exp_mask258 ) >> Exp_shift258 )) { Branch (1240:9): [True: 258, False: 1.19k]
Branch (1240:22): [True: 204, False: 54]
|
1241 | /* rv/2^bc->scale is subnormal */ |
1242 | word0(&u) = (P+2)*Exp_msk1; |
1243 | word1(&u) = 0; |
1244 | return u.d; |
1245 | } |
1246 | else { |
1247 | assert(word0(x) || word1(x)); /* x != 0.0 */ |
1248 | 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 | 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 | nd = bc->nd; |
1305 | nd0 = bc->nd0; |
1306 | p5 = nd + bc->e0; |
1307 | b = sd2b(rv, bc->scale, &p2); |
1308 | if (b == NULL) Branch (1308:9): [True: 0, False: 3.27k]
|
1309 | 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 | 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 | b = lshift(b, 1); |
1318 | if (b == NULL) Branch (1318:9): [True: 0, False: 3.27k]
|
1319 | return -1; |
1320 | b->x[0] |= 1; |
1321 | p2--; |
1322 | |
1323 | p2 -= p5; |
1324 | d = i2b(1); |
1325 | if (d == NULL) { Branch (1325:9): [True: 0, False: 3.27k]
|
1326 | Bfree(b); |
1327 | return -1; |
1328 | } |
1329 | /* Arrange for convenient computation of quotients: |
1330 | * shift left if necessary so divisor has 4 leading 0 bits. |
1331 | */ |
1332 | if (p5 > 0) { Branch (1332:9): [True: 1.07k, False: 2.19k]
|
1333 | d = pow5mult(d, p5); |
1334 | if (d == NULL) { Branch (1334:13): [True: 0, False: 1.07k]
|
1335 | Bfree(b); |
1336 | return -1; |
1337 | } |
1338 | } |
1339 | else if (p5 < 0) { Branch (1339:14): [True: 1.85k, False: 341]
|
1340 | b = pow5mult(b, -p5); |
1341 | if (b == NULL) { Branch (1341:13): [True: 0, False: 1.85k]
|
1342 | Bfree(d); |
1343 | return -1; |
1344 | } |
1345 | } |
1346 | if (p2 > 0) { Branch (1346:9): [True: 719, False: 2.55k]
|
1347 | b2 = p2; |
1348 | d2 = 0; |
1349 | } |
1350 | else { |
1351 | b2 = 0; |
1352 | d2 = -p2; |
1353 | } |
1354 | i = dshift(d, d2); |
1355 | if ((b2 += i) > 0) { Branch (1355:9): [True: 3.25k, False: 18]
|
1356 | b = lshift(b, b2); |
1357 | if (b == NULL) { Branch (1357:13): [True: 0, False: 3.25k]
|
1358 | Bfree(d); |
1359 | return -1; |
1360 | } |
1361 | } |
1362 | if ((d2 += i) > 0) { Branch (1362:9): [True: 3.25k, False: 19]
|
1363 | d = lshift(d, d2); |
1364 | if (d == NULL) { Branch (1364:13): [True: 0, False: 3.25k]
|
1365 | Bfree(b); |
1366 | 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 | if (cmp(b, d) >= 0) Branch (1373:9): [True: 72, False: 3.19k]
|
1374 | /* b/d >= 1 */ |
1375 | dd = -1; |
1376 | else { |
1377 | i = 0; |
1378 | for(;;) { |
1379 | b = multadd(b, 10, 0); |
1380 | if (b == NULL) { Branch (1380:17): [True: 0, False: 322k]
|
1381 | Bfree(d); |
1382 | return -1; |
1383 | } |
1384 | dd = s0[i < nd0 ? i321k : i+1186 ] - '0' - quorem(b, d); Branch (1384:21): [True: 321k, False: 186]
|
1385 | i++; |
1386 | |
1387 | if (dd) Branch (1387:17): [True: 2.25k, False: 319k]
|
1388 | break; |
1389 | if (!b->x[0] && b->wds == 1210k ) { Branch (1389:17): [True: 210k, False: 109k]
Branch (1389:29): [True: 945, False: 209k]
|
1390 | /* b/d == 0 */ |
1391 | dd = i < nd; |
1392 | break; |
1393 | } |
1394 | if (!(i < nd)) { Branch (1394:17): [True: 1, False: 318k]
|
1395 | /* b/d != 0, but digits of s0 exhausted */ |
1396 | dd = -1; |
1397 | break; |
1398 | } |
1399 | } |
1400 | } |
1401 | Bfree(b); |
1402 | Bfree(d); |
1403 | if (dd > 0 || (2.93k dd == 02.93k && odd943 )) Branch (1403:9): [True: 339, False: 2.93k]
Branch (1403:20): [True: 943, False: 1.98k]
Branch (1403:31): [True: 495, False: 448]
|
1404 | dval(rv) += sulp(rv, bc); |
1405 | 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 | _Py_dg_stdnan(int sign) |
1417 | { |
1418 | U rv; |
1419 | word0(&rv) = NAN_WORD0; |
1420 | word1(&rv) = NAN_WORD1; |
1421 | if (sign) Branch (1421:9): [True: 19, False: 1.97k]
|
1422 | word0(&rv) |= Sign_bit; |
1423 | 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 | _Py_dg_infinity(int sign) |
1431 | { |
1432 | U rv; |
1433 | word0(&rv) = POSINF_WORD0; |
1434 | word1(&rv) = POSINF_WORD1; |
1435 | return sign ? -1.22k dval1.22k (&rv) : dval1.71k (&rv); Branch (1435:12): [True: 1.22k, False: 1.71k]
|
1436 | } |
1437 | |
1438 | double |
1439 | _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 | Bigint *bb = NULL, *bd = NULL, *bd0 = NULL, *bs = NULL, *delta = NULL; |
1450 | size_t ndigits, fraclen; |
1451 | double result; |
1452 | |
1453 | dval(&rv) = 0.; |
1454 | |
1455 | /* Start parsing. */ |
1456 | c = *(s = s00); |
1457 | |
1458 | /* Parse optional sign, if present. */ |
1459 | sign = 0; |
1460 | switch (c) { Branch (1460:13): [True: 1.26M, False: 17.2k]
|
1461 | case '-': Branch (1461:5): [True: 13.6k, False: 1.26M]
|
1462 | sign = 1; |
1463 | /* fall through */ |
1464 | case '+': Branch (1464:5): [True: 3.63k, False: 1.27M]
|
1465 | c = *++s; |
1466 | } |
1467 | |
1468 | /* Skip leading zeros: lz is true iff there were leading zeros. */ |
1469 | s1 = s; |
1470 | while (c == '0') Branch (1470:12): [True: 1.23M, False: 1.28M]
|
1471 | c = *++s; |
1472 | 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 | s0 = s1 = s; |
1479 | while ('0' <= c && c <= '9'3.97M ) Branch (1479:12): [True: 3.97M, False: 40.1k]
Branch (1479:24): [True: 2.72M, False: 1.24M]
|
1480 | c = *++s; |
1481 | ndigits = s - s1; |
1482 | fraclen = 0; |
1483 | |
1484 | /* Parse decimal point and following digits. */ |
1485 | if (c == '.') { Branch (1485:9): [True: 35.3k, False: 1.24M]
|
1486 | c = *++s; |
1487 | if (!ndigits) { Branch (1487:13): [True: 14.0k, False: 21.2k]
|
1488 | s1 = s; |
1489 | while (c == '0') Branch (1489:20): [True: 68.0k, False: 14.0k]
|
1490 | c = *++s; |
1491 | lz = lz || s != s11.11k ; Branch (1491:18): [True: 12.8k, False: 1.11k]
Branch (1491:24): [True: 346, False: 766]
|
1492 | fraclen += (s - s1); |
1493 | s0 = s; |
1494 | } |
1495 | s1 = s; |
1496 | while ('0' <= c && c <= '9'270k ) Branch (1496:16): [True: 270k, False: 27.7k]
Branch (1496:28): [True: 262k, False: 7.58k]
|
1497 | c = *++s; |
1498 | ndigits += s - s1; |
1499 | 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 | if (!ndigits && !lz1.22M ) { Branch (1505:9): [True: 1.22M, False: 56.1k]
Branch (1505:21): [True: 5.99k, False: 1.21M]
|
1506 | if (se) Branch (1506:13): [True: 5.99k, False: 0]
|
1507 | *se = (char *)s00; |
1508 | 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 | if (ndigits > MAX_DIGITS || fraclen > MAX_DIGITS) { Branch (1513:9): [True: 0, False: 1.27M]
Branch (1513:33): [True: 0, False: 1.27M]
|
1514 | if (se) Branch (1514:13): [True: 0, False: 0]
|
1515 | *se = (char *)s00; |
1516 | goto parse_error; |
1517 | } |
1518 | nd = (int)ndigits; |
1519 | nd0 = (int)ndigits - (int)fraclen; |
1520 | |
1521 | /* Parse exponent. */ |
1522 | e = 0; |
1523 | if (c == 'e' || c == 'E'37.3k ) { Branch (1523:9): [True: 1.23M, False: 37.3k]
Branch (1523:21): [True: 4.30k, False: 33.0k]
|
1524 | s00 = s; |
1525 | c = *++s; |
1526 | |
1527 | /* Exponent sign. */ |
1528 | esign = 0; |
1529 | switch (c) { Branch (1529:17): [True: 8.45k, False: 1.23M]
|
1530 | case '-': Branch (1530:9): [True: 1.23M, False: 11.9k]
|
1531 | esign = 1; |
1532 | /* fall through */ |
1533 | case '+': Branch (1533:9): [True: 3.51k, False: 1.23M]
|
1534 | c = *++s; |
1535 | } |
1536 | |
1537 | /* Skip zeros. lz is true iff there are leading zeros. */ |
1538 | s1 = s; |
1539 | while (c == '0') Branch (1539:16): [True: 5.01k, False: 1.24M]
|
1540 | c = *++s; |
1541 | lz = s != s1; |
1542 | |
1543 | /* Get absolute value of the exponent. */ |
1544 | s1 = s; |
1545 | abs_exp = 0; |
1546 | while ('0' <= c && c <= '9'1.27M ) { Branch (1546:16): [True: 1.27M, False: 1.24M]
Branch (1546:28): [True: 1.27M, False: 438]
|
1547 | abs_exp = 10*abs_exp + (c - '0'); |
1548 | 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 | if (s - s1 > 9 || abs_exp > MAX_ABS_EXP) Branch (1554:13): [True: 0, False: 1.24M]
Branch (1554:27): [True: 0, False: 1.24M]
|
1555 | e = (int)MAX_ABS_EXP; |
1556 | else |
1557 | e = (int)abs_exp; |
1558 | if (esign) Branch (1558:13): [True: 1.23M, False: 11.9k]
|
1559 | e = -e; |
1560 | |
1561 | /* A valid exponent must have at least one digit. */ |
1562 | if (s == s1 && !lz3.20k ) Branch (1562:13): [True: 3.20k, False: 1.23M]
Branch (1562:24): [True: 1, False: 3.20k]
|
1563 | s = s00; |
1564 | } |
1565 | |
1566 | /* Adjust exponent to take into account position of the point. */ |
1567 | e -= nd - nd0; |
1568 | if (nd0 <= 0) Branch (1568:9): [True: 1.22M, False: 47.7k]
|
1569 | nd0 = nd; |
1570 | |
1571 | /* Finished parsing. Set se to indicate how far we parsed */ |
1572 | if (se) Branch (1572:9): [True: 58.2k, False: 1.21M]
|
1573 | *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 | if (!nd) Branch (1577:9): [True: 1.21M, False: 56.1k]
|
1578 | goto ret; |
1579 | for (i = nd; 56.1k i > 0; ) { Branch (1579:18): [True: 242k, False: 0]
|
1580 | --i; |
1581 | if (s0[i < nd0 ? i212k : i+129.8k ] != '0') { Branch (1581:13): [True: 56.1k, False: 186k]
Branch (1581:16): [True: 212k, False: 29.8k]
|
1582 | ++i; |
1583 | break; |
1584 | } |
1585 | } |
1586 | e += nd - i; |
1587 | nd = i; |
1588 | if (nd0 > nd) Branch (1588:9): [True: 5.98k, False: 50.1k]
|
1589 | 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 | bc.e0 = e1 = e; |
1628 | y = z = 0; |
1629 | for (i = 0; i < nd; i++444k ) { Branch (1629:17): [True: 462k, False: 38.5k]
|
1630 | if (i < 9) Branch (1630:13): [True: 296k, False: 165k]
|
1631 | y = 10*y + s0[i < nd0 ? i215k : i+181.6k ] - '0'; Branch (1631:27): [True: 215k, False: 81.6k]
|
1632 | else if (i < DBL_DIG+1) Branch (1632:18): [True: 147k, False: 17.5k]
|
1633 | z = 10*z + s0[i < nd0 ? i90.5k : i+157.3k ] - '0'; Branch (1633:27): [True: 90.5k, False: 57.3k]
|
1634 | else |
1635 | break; |
1636 | } |
1637 | |
1638 | k = nd < DBL_DIG + 1 ? nd37.1k : DBL_DIG + 118.9k ; Branch (1638:9): [True: 37.1k, False: 18.9k]
|
1639 | dval(&rv) = y; |
1640 | if (k > 9) { Branch (1640:9): [True: 23.7k, False: 32.3k]
|
1641 | dval(&rv) = tens[k - 9] * dval(&rv) + z; |
1642 | } |
1643 | if (nd <= DBL_DIG Branch (1643:9): [True: 37.1k, False: 18.9k]
|
1644 | && Flt_Rounds37.1k == 137.1k Branch (1644:12): [True: 37.1k, False: 0]
|
1645 | ) { |
1646 | if (!e) Branch (1646:13): [True: 7.16k, False: 30.0k]
|
1647 | goto ret; |
1648 | if (e > 0) { Branch (1648:13): [True: 9.14k, False: 20.8k]
|
1649 | if (e <= Ten_pmax) { Branch (1649:17): [True: 5.62k, False: 3.52k]
|
1650 | dval(&rv) *= tens[e]; |
1651 | goto ret; |
1652 | } |
1653 | i = DBL_DIG - nd; |
1654 | if (e <= Ten_pmax + i) { Branch (1654:17): [True: 799, False: 2.72k]
|
1655 | /* A fancier test would sometimes let us do |
1656 | * this for larger i values. |
1657 | */ |
1658 | e -= i; |
1659 | dval(&rv) *= tens[i]; |
1660 | dval(&rv) *= tens[e]; |
1661 | goto ret; |
1662 | } |
1663 | } |
1664 | else if (e >= -Ten_pmax) { Branch (1664:18): [True: 17.8k, False: 3.00k]
|
1665 | dval(&rv) /= tens[-e]; |
1666 | goto ret; |
1667 | } |
1668 | } |
1669 | e1 += nd - k; |
1670 | |
1671 | bc.scale = 0; |
1672 | |
1673 | /* Get starting approximation = rv * 10**e1 */ |
1674 | |
1675 | if (e1 > 0) { Branch (1675:9): [True: 7.40k, False: 17.2k]
|
1676 | if ((i = e1 & 15)) Branch (1676:13): [True: 7.13k, False: 270]
|
1677 | dval(&rv) *= tens[i]; |
1678 | if (e1 &= ~15) { Branch (1678:13): [True: 6.19k, False: 1.21k]
|
1679 | if (e1 > DBL_MAX_10_EXP) Branch (1679:17): [True: 687, False: 5.50k]
|
1680 | goto ovfl; |
1681 | e1 >>= 4; |
1682 | for(j = 0; e1 > 1; j++, e1 >>= 113.3k ) Branch (1682:24): [True: 13.3k, False: 5.50k]
|
1683 | if (e1 & 1) Branch (1683:21): [True: 4.65k, False: 8.66k]
|
1684 | dval(&rv) *= bigtens[j]; |
1685 | /* The last multiplication could overflow. */ |
1686 | word0(&rv) -= P*Exp_msk1; |
1687 | dval(&rv) *= bigtens[j]; |
1688 | if ((z = word0(&rv) & Exp_mask) Branch (1688:17): [True: 64, False: 5.44k]
|
1689 | > Exp_msk1*(DBL_MAX_EXP+Bias-P)) |
1690 | goto ovfl; |
1691 | if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) { Branch (1691:17): [True: 417, False: 5.02k]
|
1692 | /* set to largest number */ |
1693 | /* (Can't trust DBL_MAX) */ |
1694 | word0(&rv) = Big0; |
1695 | word1(&rv) = Big1; |
1696 | } |
1697 | else |
1698 | word0(&rv) += P*Exp_msk1; |
1699 | } |
1700 | } |
1701 | else if (e1 < 0) { Branch (1701:14): [True: 16.9k, False: 354]
|
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 | e1 = -e1; |
1713 | if ((i = e1 & 15)) Branch (1713:13): [True: 14.4k, False: 2.41k]
|
1714 | dval(&rv) /= tens[i]; |
1715 | if (e1 >>= 4) { Branch (1715:13): [True: 10.8k, False: 6.03k]
|
1716 | if (e1 >= 1 << n_bigtens) Branch (1716:17): [True: 118, False: 10.7k]
|
1717 | goto undfl; |
1718 | if (e1 & Scale_Bit) Branch (1718:17): [True: 4.46k, False: 6.30k]
|
1719 | bc.scale = 2*P; |
1720 | for(j = 0; e1 > 0; j++, e1 >>= 135.6k ) Branch (1720:24): [True: 35.6k, False: 10.7k]
|
1721 | if (e1 & 1) Branch (1721:21): [True: 20.7k, False: 14.8k]
|
1722 | dval(&rv) *= tinytens[j]; |
1723 | if (bc.scale && (j = 2*4.46k P4.46k + 1 - ((word04.46k (&rv) & Exp_mask4.46k ) Branch (1723:17): [True: 4.46k, False: 6.30k]
Branch (1723:29): [True: 3.49k, False: 969]
|
1724 | >> Exp_shift)) > 0) { |
1725 | /* scaled rv is denormal; clear j low bits */ |
1726 | if (j >= 32) { Branch (1726:21): [True: 2.25k, False: 1.23k]
|
1727 | word1(&rv) = 0; |
1728 | if (j >= 53) Branch (1728:25): [True: 1.23k, False: 1.02k]
|
1729 | word0(&rv) = (P+2)*Exp_msk1; |
1730 | else |
1731 | word0(&rv) &= 0xffffffff << (j-32); |
1732 | } |
1733 | else |
1734 | word1(&rv) &= 0xffffffff << j; |
1735 | } |
1736 | if (!dval(&rv)) Branch (1736:17): [True: 0, False: 10.7k]
|
1737 | 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 | bc.nd = nd; |
1746 | 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 | if (nd > STRTOD_DIGLIM) { Branch (1749:9): [True: 6.63k, False: 17.1k]
|
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 | for (i = 18; i > 0; ) { Branch (1757:22): [True: 7.52k, False: 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 | --i; |
1761 | if (s0[i < nd0 ? i6.22k : i+11.30k ] != '0') { Branch (1761:17): [True: 6.63k, False: 893]
Branch (1761:20): [True: 6.22k, False: 1.30k]
|
1762 | ++i; |
1763 | break; |
1764 | } |
1765 | } |
1766 | e += nd - i; |
1767 | nd = i; |
1768 | if (nd0 > nd) Branch (1768:13): [True: 5.46k, False: 1.16k]
|
1769 | nd0 = nd; |
1770 | if (nd < 9) { /* must recompute y */ Branch (1770:13): [True: 14, False: 6.61k]
|
1771 | y = 0; |
1772 | for(i = 0; i < nd0; ++i15 ) Branch (1772:24): [True: 15, False: 14]
|
1773 | y = 10*y + s0[i] - '0'; |
1774 | for(; i < nd; ++i1 ) Branch (1774:19): [True: 1, False: 14]
|
1775 | y = 10*y + s0[i+1] - '0'; |
1776 | } |
1777 | } |
1778 | bd0 = s2b(s0, nd0, nd, y); |
1779 | if (bd0 == NULL) Branch (1779:9): [True: 0, False: 23.8k]
|
1780 | 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(;;)23.8k { |
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 | bd = Balloc(bd0->k); |
1811 | if (bd == NULL) { Branch (1811:13): [True: 0, False: 30.4k]
|
1812 | goto failed_malloc; |
1813 | } |
1814 | Bcopy(bd, bd0); |
1815 | bb = sd2b(&rv, bc.scale, &bbe); /* srv = bb * 2^bbe */ |
1816 | if (bb == NULL) { Branch (1816:13): [True: 0, False: 30.4k]
|
1817 | 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 | odd = bb->x[0] & 1; |
1822 | |
1823 | /* tdv = bd * 10**e; srv = bb * 2**bbe */ |
1824 | bs = i2b(1); |
1825 | if (bs == NULL) { Branch (1825:13): [True: 0, False: 30.4k]
|
1826 | goto failed_malloc; |
1827 | } |
1828 | |
1829 | if (e >= 0) { Branch (1829:13): [True: 8.39k, False: 22.0k]
|
1830 | bb2 = bb5 = 0; |
1831 | bd2 = bd5 = e; |
1832 | } |
1833 | else { |
1834 | bb2 = bb5 = -e; |
1835 | bd2 = bd5 = 0; |
1836 | } |
1837 | if (bbe >= 0) Branch (1837:13): [True: 8.48k, False: 21.9k]
|
1838 | bb2 += bbe; |
1839 | else |
1840 | bd2 -= bbe; |
1841 | bs2 = bb2; |
1842 | bb2++; |
1843 | 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 | i = bb2 < bd2 ? bb221.3k : bd29.09k ; Branch (1863:13): [True: 21.3k, False: 9.09k]
|
1864 | if (i > bs2) Branch (1864:13): [True: 21.3k, False: 9.09k]
|
1865 | i = bs2; |
1866 | if (i > 0) { Branch (1866:13): [True: 30.4k, False: 12]
|
1867 | bb2 -= i; |
1868 | bd2 -= i; |
1869 | bs2 -= i; |
1870 | } |
1871 | |
1872 | /* Scale bb, bd, bs by the appropriate powers of 2 and 5. */ |
1873 | if (bb5 > 0) { Branch (1873:13): [True: 22.0k, False: 8.39k]
|
1874 | bs = pow5mult(bs, bb5); |
1875 | if (bs == NULL) { Branch (1875:17): [True: 0, False: 22.0k]
|
1876 | goto failed_malloc; |
1877 | } |
1878 | Bigint *bb1 = mult(bs, bb); |
1879 | Bfree(bb); |
1880 | bb = bb1; |
1881 | if (bb == NULL) { Branch (1881:17): [True: 0, False: 22.0k]
|
1882 | goto failed_malloc; |
1883 | } |
1884 | } |
1885 | if (bb2 > 0) { Branch (1885:13): [True: 30.4k, False: 0]
|
1886 | bb = lshift(bb, bb2); |
1887 | if (bb == NULL) { Branch (1887:17): [True: 0, False: 30.4k]
|
1888 | goto failed_malloc; |
1889 | } |
1890 | } |
1891 | if (bd5 > 0) { Branch (1891:13): [True: 7.83k, False: 22.6k]
|
1892 | bd = pow5mult(bd, bd5); |
1893 | if (bd == NULL) { Branch (1893:17): [True: 0, False: 7.83k]
|
1894 | goto failed_malloc; |
1895 | } |
1896 | } |
1897 | if (bd2 > 0) { Branch (1897:13): [True: 21.3k, False: 9.09k]
|
1898 | bd = lshift(bd, bd2); |
1899 | if (bd == NULL) { Branch (1899:17): [True: 0, False: 21.3k]
|
1900 | goto failed_malloc; |
1901 | } |
1902 | } |
1903 | if (bs2 > 0) { Branch (1903:13): [True: 8.07k, False: 22.3k]
|
1904 | bs = lshift(bs, bs2); |
1905 | if (bs == NULL) { Branch (1905:17): [True: 0, False: 8.07k]
|
1906 | 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 | delta = diff(bb, bd); |
1915 | if (delta == NULL) { Branch (1915:13): [True: 0, False: 30.4k]
|
1916 | goto failed_malloc; |
1917 | } |
1918 | dsign = delta->sign; |
1919 | delta->sign = 0; |
1920 | i = cmp(delta, bs); |
1921 | if (bc.nd > nd && i <= 011.1k ) { Branch (1921:13): [True: 11.1k, False: 19.3k]
Branch (1921:27): [True: 6.06k, False: 5.05k]
|
1922 | if (dsign) Branch (1922:17): [True: 3.10k, False: 2.95k]
|
1923 | 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 | if (!word1(&rv) && !(756 word0756 (&rv) & Bndry_mask756 )) { Branch (1936:17): [True: 756, False: 2.19k]
Branch (1936:32): [True: 577, False: 179]
|
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 | 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 | if (j - bc.scale >= 2) { Branch (1942:21): [True: 162, False: 415]
|
1943 | dval(&rv) -= 0.5 * sulp(&rv, &bc); |
1944 | break; /* Use bigcomp. */ |
1945 | } |
1946 | } |
1947 | |
1948 | { |
1949 | bc.nd = nd; |
1950 | i = -1; /* Discarded digits make delta smaller. */ |
1951 | } |
1952 | } |
1953 | |
1954 | if (i < 0) { Branch (1954:13): [True: 10.1k, False: 17.0k]
|
1955 | /* Error is less than half an ulp -- check for |
1956 | * special case of mantissa a power of two. |
1957 | */ |
1958 | if (dsign || word16.49k (&rv) || word0957 (&rv) & 957 Bndry_mask957 Branch (1958:17): [True: 3.68k, False: 6.49k]
Branch (1958:40): [True: 288, False: 669]
|
1959 | || (669 word0669 (&rv) & Exp_mask669 ) <= (2*P669 +1)*Exp_msk1669 Branch (1959:20): [True: 631, False: 38]
|
1960 | ) { |
1961 | break; |
1962 | } |
1963 | if (!delta->x[0] && delta->wds <= 137 ) { Branch (1963:17): [True: 37, False: 1]
Branch (1963:33): [True: 2, False: 35]
|
1964 | /* exact result */ |
1965 | break; |
1966 | } |
1967 | delta = lshift(delta,Log2P); |
1968 | if (delta == NULL) { Branch (1968:17): [True: 0, False: 36]
|
1969 | goto failed_malloc; |
1970 | } |
1971 | if (cmp(delta, bs) > 0) Branch (1971:17): [True: 3, False: 33]
|
1972 | goto drop_down; |
1973 | break; |
1974 | } |
1975 | if (i == 0) { Branch (1975:13): [True: 2.34k, False: 14.6k]
|
1976 | /* exactly half-way between */ |
1977 | if (dsign) { Branch (1977:17): [True: 1.62k, False: 720]
|
1978 | if ((word0(&rv) & Bndry_mask1) == Bndry_mask1 Branch (1978:21): [True: 0, False: 1.62k]
|
1979 | && word10 (&rv) == ( Branch (1979:25): [True: 0, False: 0]
|
1980 | (bc.scale && Branch (1980:26): [True: 0, False: 0]
|
1981 | (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1) ? Branch (1981:26): [True: 0, False: 0]
|
1982 | (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) : |
1983 | 0xffffffff)) { |
1984 | /*boundary case -- increment exponent*/ |
1985 | word0(&rv) = (word0(&rv) & Exp_mask) |
1986 | + Exp_msk1 |
1987 | ; |
1988 | word1(&rv) = 0; |
1989 | /* dsign = 0; */ |
1990 | break; |
1991 | } |
1992 | } |
1993 | else if (!(word0(&rv) & Bndry_mask) && !0 word10 (&rv)) { Branch (1993:22): [True: 0, False: 720]
Branch (1993:52): [True: 0, False: 0]
|
1994 | drop_down: |
1995 | /* boundary case -- decrement exponent */ |
1996 | if (bc.scale) { Branch (1996:21): [True: 0, False: 3]
|
1997 | L = word0(&rv) & Exp_mask; |
1998 | if (L <= (2*P+1)*Exp_msk1) { Branch (1998:25): [True: 0, False: 0]
|
1999 | if (L > (P+2)*Exp_msk1) Branch (1999:29): [True: 0, False: 0]
|
2000 | /* round even ==> */ |
2001 | /* accept rv */ |
2002 | break; |
2003 | /* rv = smallest denormal */ |
2004 | if (bc.nd > nd) Branch (2004:29): [True: 0, False: 0]
|
2005 | break; |
2006 | goto undfl; |
2007 | } |
2008 | } |
2009 | L = (word0(&rv) & Exp_mask) - Exp_msk1; |
2010 | word0(&rv) = L | Bndry_mask1; |
2011 | word1(&rv) = 0xffffffff; |
2012 | break; |
2013 | } |
2014 | if (!odd) Branch (2014:17): [True: 1.89k, False: 453]
|
2015 | break; |
2016 | if (dsign) Branch (2016:17): [True: 447, False: 6]
|
2017 | dval(&rv) += sulp(&rv, &bc); |
2018 | else { |
2019 | dval(&rv) -= sulp(&rv, &bc); |
2020 | if (!dval(&rv)) { Branch (2020:21): [True: 0, False: 6]
|
2021 | if (bc.nd >nd) Branch (2021:25): [True: 0, False: 0]
|
2022 | break; |
2023 | goto undfl; |
2024 | } |
2025 | } |
2026 | /* dsign = 1 - dsign; */ |
2027 | break; |
2028 | } |
2029 | if ((aadj = ratio(delta, bs)) <= 2.) { Branch (2029:13): [True: 6.20k, False: 8.44k]
|
2030 | if (dsign) Branch (2030:17): [True: 3.86k, False: 2.33k]
|
2031 | aadj = aadj1 = 1.; |
2032 | else if (word1(&rv) || word0862 (&rv) & 862 Bndry_mask862 ) { Branch (2032:36): [True: 0, False: 862]
|
2033 | if (word1(&rv) == Tiny1 && !0 word00 (&rv)) { Branch (2033:21): [True: 0, False: 1.47k]
Branch (2033:44): [True: 0, False: 0]
|
2034 | if (bc.nd >nd) Branch (2034:25): [True: 0, False: 0]
|
2035 | break; |
2036 | goto undfl; |
2037 | } |
2038 | aadj = 1.; |
2039 | aadj1 = -1.; |
2040 | } |
2041 | else { |
2042 | /* special case -- power of FLT_RADIX to be */ |
2043 | /* rounded down... */ |
2044 | |
2045 | if (aadj < 2./FLT_RADIX) Branch (2045:21): [True: 0, False: 862]
|
2046 | aadj = 1./FLT_RADIX; |
2047 | else |
2048 | aadj *= 0.5; |
2049 | aadj1 = -aadj; |
2050 | } |
2051 | } |
2052 | else { |
2053 | aadj *= 0.5; |
2054 | aadj1 = dsign ? aadj7.91k : -aadj534 ; Branch (2054:21): [True: 7.91k, False: 534]
|
2055 | if (Flt_Rounds == 0) Branch (2055:17): [True: 0, False: 8.44k]
|
2056 | aadj1 += 0.5; |
2057 | } |
2058 | y = word0(&rv) & Exp_mask; |
2059 | |
2060 | /* Check for overflow */ |
2061 | |
2062 | if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) { Branch (2062:13): [True: 1.42k, False: 13.2k]
|
2063 | dval(&rv0) = dval(&rv); |
2064 | word0(&rv) -= P*Exp_msk1; |
2065 | adj.d = aadj1 * ulp(&rv); |
2066 | dval(&rv) += adj.d; |
2067 | if ((word0(&rv) & Exp_mask) >= Branch (2067:17): [True: 761, False: 664]
|
2068 | Exp_msk1*(DBL_MAX_EXP+Bias-P)) { |
2069 | if (word0(&rv0) == Big0 && word1(&rv0) == Big1) { Branch (2069:21): [True: 761, False: 0]
Branch (2069:44): [True: 589, False: 172]
|
2070 | goto ovfl; |
2071 | } |
2072 | word0(&rv) = Big0; |
2073 | word1(&rv) = Big1; |
2074 | goto cont; |
2075 | } |
2076 | else |
2077 | word0(&rv) += P*Exp_msk1; |
2078 | } |
2079 | else { |
2080 | if (bc.scale && y <= 2*2.96k P2.96k *Exp_msk12.96k ) { Branch (2080:17): [True: 2.96k, False: 10.2k]
Branch (2080:29): [True: 2.37k, False: 588]
|
2081 | if (aadj <= 0x7fffffff) { Branch (2081:21): [True: 2.37k, False: 0]
|
2082 | if ((z = (ULong)aadj) <= 0) Branch (2082:25): [True: 748, False: 1.62k]
|
2083 | z = 1; |
2084 | aadj = z; |
2085 | aadj1 = dsign ? aadj1.51k : -aadj862 ; Branch (2085:29): [True: 1.51k, False: 862]
|
2086 | } |
2087 | dval(&aadj2) = aadj1; |
2088 | word0(&aadj2) += (2*P+1)*Exp_msk1 - y; |
2089 | aadj1 = dval(&aadj2); |
2090 | } |
2091 | adj.d = aadj1 * ulp(&rv); |
2092 | dval(&rv) += adj.d; |
2093 | } |
2094 | z = word0(&rv) & Exp_mask; |
2095 | if (bc.nd == nd) { Branch (2095:13): [True: 9.57k, False: 4.32k]
|
2096 | if (!bc.scale) Branch (2096:17): [True: 8.45k, False: 1.12k]
|
2097 | if (y == z) { Branch (2097:21): [True: 8.41k, False: 31]
|
2098 | /* Can we stop now? */ |
2099 | L = (Long)aadj; |
2100 | aadj -= L; |
2101 | /* The tolerances below are conservative. */ |
2102 | if (dsign || word11.57k (&rv) || word00 (&rv) & 0 Bndry_mask0 ) { Branch (2102:25): [True: 6.84k, False: 1.57k]
Branch (2102:48): [True: 0, False: 0]
|
2103 | if (aadj < .4999999 || aadj > .50000012.18k ) Branch (2103:29): [True: 6.23k, False: 2.18k]
Branch (2103:48): [True: 1.18k, False: 1.00k]
|
2104 | break; |
2105 | } |
2106 | else if (aadj < .4999999/FLT_RADIX) Branch (2106:30): [True: 0, False: 0]
|
2107 | break; |
2108 | } |
2109 | } |
2110 | cont: |
2111 | Bfree(bb); bb = NULL; |
2112 | Bfree(bd); bd = NULL; |
2113 | Bfree(bs); bs = NULL; |
2114 | Bfree(delta); delta = NULL; |
2115 | } |
2116 | if (bc.nd > nd) { Branch (2116:9): [True: 3.27k, False: 19.9k]
|
2117 | error = bigcomp(&rv, s0, &bc); |
2118 | if (error) Branch (2118:13): [True: 0, False: 3.27k]
|
2119 | goto failed_malloc; |
2120 | } |
2121 | |
2122 | if (bc.scale) { Branch (2122:9): [True: 4.46k, False: 18.7k]
|
2123 | word0(&rv0) = Exp_1 - 2*P*Exp_msk1; |
2124 | word1(&rv0) = 0; |
2125 | dval(&rv) *= dval(&rv0); |
2126 | } |
2127 | |
2128 | ret: |
2129 | result = sign ? -11.8k dval11.8k (&rv) : dval1.26M (&rv); Branch (2129:14): [True: 11.8k, False: 1.26M]
|
2130 | goto done; |
2131 | |
2132 | parse_error: |
2133 | result = 0.0; |
2134 | goto done; |
2135 | |
2136 | failed_malloc: |
2137 | errno = ENOMEM; |
2138 | result = -1.0; |
2139 | goto done; |
2140 | |
2141 | undfl: |
2142 | result = sign ? -0.050 : 0.068 ; Branch (2142:14): [True: 50, False: 68]
|
2143 | goto done; |
2144 | |
2145 | ovfl: |
2146 | errno = ERANGE; |
2147 | /* Can't trust HUGE_VAL */ |
2148 | word0(&rv) = Exp_mask; |
2149 | word1(&rv) = 0; |
2150 | result = sign ? -123 dval123 (&rv) : dval1.21k (&rv); Branch (2150:14): [True: 123, False: 1.21k]
|
2151 | goto done; |
2152 | |
2153 | done: |
2154 | Bfree(bb); |
2155 | Bfree(bd); |
2156 | Bfree(bs); |
2157 | Bfree(bd0); |
2158 | Bfree(delta); |
2159 | return result; |
2160 | |
2161 | } |
2162 | |
2163 | static char * |
2164 | rv_alloc(int i) |
2165 | { |
2166 | int j, k, *r; |
2167 | |
2168 | j = sizeof(ULong); |
2169 | for(k = 0; |
2170 | sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= (unsigned)i; Branch (2170:9): [True: 51.2k, False: 1.41M]
|
2171 | j <<= 151.2k ) |
2172 | k++; |
2173 | r = (int*)Balloc(k); |
2174 | if (r == NULL) Branch (2174:9): [True: 0, False: 1.41M]
|
2175 | return NULL; |
2176 | *r = k; |
2177 | return (char *)(r+1); |
2178 | } |
2179 | |
2180 | static char * |
2181 | nrv_alloc(const char *s, char **rve, int n) |
2182 | { |
2183 | char *rv, *t; |
2184 | |
2185 | rv = rv_alloc(n); |
2186 | if (rv == NULL) Branch (2186:9): [True: 0, False: 20.6k]
|
2187 | return NULL; |
2188 | t = rv; |
2189 | while((*t = *s++)) t++120k ; Branch (2189:11): [True: 120k, False: 20.6k]
|
2190 | if (rve) Branch (2190:9): [True: 20.6k, False: 0]
|
2191 | *rve = t; |
2192 | 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 | _Py_dg_freedtoa(char *s) |
2203 | { |
2204 | Bigint *b = (Bigint *)((int *)s - 1); |
2205 | b->maxwds = 1 << (b->k = *(int*)b); |
2206 | Bfree(b); |
2207 | } |
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 | _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 | mlo = mhi = S = 0; |
2299 | s0 = 0; |
2300 | |
2301 | u.d = dd; |
2302 | if (word0(&u) & Sign_bit) { Branch (2302:9): [True: 75.8k, False: 1.34M]
|
2303 | /* set sign for everything, including 0's and NaNs */ |
2304 | *sign = 1; |
2305 | word0(&u) &= ~Sign_bit; /* clear sign bit */ |
2306 | } |
2307 | else |
2308 | *sign = 0; |
2309 | |
2310 | /* quick return for Infinities, NaNs and zeros */ |
2311 | if ((word0(&u) & Exp_mask) == Exp_mask) Branch (2311:9): [True: 15.3k, False: 1.40M]
|
2312 | { |
2313 | /* Infinity or NaN */ |
2314 | *decpt = 9999; |
2315 | if (!word1(&u) && !(word0(&u) & 0xfffff)) Branch (2315:13): [True: 15.3k, False: 0]
Branch (2315:27): [True: 13.8k, False: 1.53k]
|
2316 | return nrv_alloc("Infinity", rve, 8); |
2317 | return nrv_alloc("NaN", rve, 3); |
2318 | } |
2319 | if (!dval(&u)) { Branch (2319:9): [True: 5.22k, False: 1.39M]
|
2320 | *decpt = 1; |
2321 | 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 | b = d2b(&u, &be, &bbits); |
2327 | if (b == NULL) Branch (2327:9): [True: 0, False: 1.39M]
|
2328 | goto failed_malloc; |
2329 | if ((i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1)))) { Branch (2329:9): [True: 1.39M, False: 991]
|
2330 | dval(&d2) = dval(&u); |
2331 | word0(&d2) &= Frac_mask1; |
2332 | 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 | i -= Bias; |
2357 | denorm = 0; |
2358 | } |
2359 | else { |
2360 | /* d is denormalized */ |
2361 | |
2362 | i = bbits + be + (Bias + (P-1) - 1); |
2363 | x = i > 32 ? word0343 (&u) << (64 - i) | 343 word1343 (&u) >> (i - 32) Branch (2363:13): [True: 343, False: 648]
|
2364 | : word1648 (&u) << (32 - i)648 ; |
2365 | dval(&d2) = x; |
2366 | word0(&d2) -= 31*Exp_msk1; /* adjust exponent */ |
2367 | i -= (Bias + (P-1) - 1) + 1; |
2368 | denorm = 1; |
2369 | } |
2370 | ds = (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 + |
2371 | i*0.301029995663981; |
2372 | k = (int)ds; |
2373 | if (ds < 0. && ds != k1.29M ) Branch (2373:9): [True: 1.29M, False: 109k]
Branch (2373:20): [True: 1.29M, False: 0]
|
2374 | k--; /* want k = floor(ds) */ |
2375 | k_check = 1; |
2376 | if (k >= 0 && k <= 109k Ten_pmax109k ) { Branch (2376:9): [True: 109k, False: 1.29M]
Branch (2376:19): [True: 64.6k, False: 44.6k]
|
2377 | if (dval(&u) < tens[k]) Branch (2377:13): [True: 1.29k, False: 63.3k]
|
2378 | k--; |
2379 | k_check = 0; |
2380 | } |
2381 | j = bbits - i - 1; |
2382 | if (j >= 0) { Branch (2382:9): [True: 1.33M, False: 63.4k]
|
2383 | b2 = 0; |
2384 | s2 = j; |
2385 | } |
2386 | else { |
2387 | b2 = -j; |
2388 | s2 = 0; |
2389 | } |
2390 | if (k >= 0) { Branch (2390:9): [True: 108k, False: 1.29M]
|
2391 | b5 = 0; |
2392 | s5 = k; |
2393 | s2 += k; |
2394 | } |
2395 | else { |
2396 | b2 -= k; |
2397 | b5 = -k; |
2398 | s5 = 0; |
2399 | } |
2400 | if (mode < 0 || mode > 9) Branch (2400:9): [True: 0, False: 1.39M]
Branch (2400:21): [True: 0, False: 1.39M]
|
2401 | mode = 0; |
2402 | |
2403 | try_quick = 1; |
2404 | |
2405 | if (mode > 5) { Branch (2405:9): [True: 0, False: 1.39M]
|
2406 | mode -= 4; |
2407 | try_quick = 0; |
2408 | } |
2409 | leftright = 1; |
2410 | ilim = ilim1 = -1; /* Values for cases 0 and 1; done here to */ |
2411 | /* silence erroneous "gcc -Wall" warning. */ |
2412 | switch(mode) { Branch (2412:12): [True: 0, False: 1.39M]
|
2413 | case 0: Branch (2413:5): [True: 144k, False: 1.25M]
|
2414 | case 1: Branch (2414:5): [True: 0, False: 1.39M]
|
2415 | i = 18; |
2416 | ndigits = 0; |
2417 | break; |
2418 | case 2: Branch (2418:5): [True: 8.77k, False: 1.39M]
|
2419 | leftright = 0; |
2420 | /* fall through */ |
2421 | case 4: Branch (2421:5): [True: 0, False: 1.39M]
|
2422 | if (ndigits <= 0) Branch (2422:13): [True: 0, False: 8.77k]
|
2423 | ndigits = 1; |
2424 | ilim = ilim1 = i = ndigits; |
2425 | break; |
2426 | case 3: Branch (2426:5): [True: 1.24M, False: 152k]
|
2427 | leftright = 0; |
2428 | /* fall through */ |
2429 | case 5: Branch (2429:5): [True: 0, False: 1.39M]
|
2430 | i = ndigits + k + 1; |
2431 | ilim = i; |
2432 | ilim1 = i - 1; |
2433 | if (i <= 0) Branch (2433:13): [True: 1.21M, False: 32.4k]
|
2434 | i = 1; |
2435 | } |
2436 | s0 = rv_alloc(i); |
2437 | if (s0 == NULL) Branch (2437:9): [True: 0, False: 1.39M]
|
2438 | goto failed_malloc; |
2439 | s = s0; |
2440 | |
2441 | |
2442 | if (ilim >= 0 && ilim <= 43.9k Quick_max43.9k && try_quick30.0k ) { Branch (2442:9): [True: 43.9k, False: 1.35M]
Branch (2442:22): [True: 30.0k, False: 13.8k]
Branch (2442:43): [True: 30.0k, False: 0]
|
2443 | |
2444 | /* Try to get by with floating-point arithmetic. */ |
2445 | |
2446 | i = 0; |
2447 | dval(&d2) = dval(&u); |
2448 | k0 = k; |
2449 | ilim0 = ilim; |
2450 | ieps = 2; /* conservative */ |
2451 | if (k > 0) { Branch (2451:13): [True: 10.3k, False: 19.7k]
|
2452 | ds = tens[k&0xf]; |
2453 | j = k >> 4; |
2454 | if (j & Bletch) { Branch (2454:17): [True: 2, False: 10.3k]
|
2455 | /* prevent overflows */ |
2456 | j &= Bletch - 1; |
2457 | dval(&u) /= bigtens[n_bigtens-1]; |
2458 | ieps++; |
2459 | } |
2460 | for(; j; j >>= 1, i++440 ) Branch (2460:19): [True: 440, False: 10.3k]
|
2461 | if (j & 1) { Branch (2461:21): [True: 284, False: 156]
|
2462 | ieps++; |
2463 | ds *= bigtens[i]; |
2464 | } |
2465 | dval(&u) /= ds; |
2466 | } |
2467 | else if ((j1 = -k)) { Branch (2467:18): [True: 14.7k, False: 4.94k]
|
2468 | dval(&u) *= tens[j1 & 0xf]; |
2469 | for(j = j1 >> 4; j; j >>= 1, i++394 ) Branch (2469:30): [True: 394, False: 14.7k]
|
2470 | if (j & 1) { Branch (2470:21): [True: 252, False: 142]
|
2471 | ieps++; |
2472 | dval(&u) *= bigtens[i]; |
2473 | } |
2474 | } |
2475 | if (k_check && dval14.5k (&u) < 1.14.5k && ilim > 07 ) { Branch (2475:13): [True: 14.5k, False: 15.5k]
Branch (2475:24): [True: 7, False: 14.5k]
Branch (2475:41): [True: 3, False: 4]
|
2476 | if (ilim1 <= 0) Branch (2476:17): [True: 1, False: 2]
|
2477 | goto fast_failed; |
2478 | ilim = ilim1; |
2479 | k--; |
2480 | dval(&u) *= 10.; |
2481 | ieps++; |
2482 | } |
2483 | dval(&eps) = ieps*dval(&u) + 7.; |
2484 | word0(&eps) -= (P-1)*Exp_msk1; |
2485 | if (ilim == 0) { Branch (2485:13): [True: 2.65k, False: 27.4k]
|
2486 | S = mhi = 0; |
2487 | dval(&u) -= 5.; |
2488 | if (dval(&u) > dval(&eps)) Branch (2488:17): [True: 803, False: 1.84k]
|
2489 | goto one_digit; |
2490 | if (dval(&u) < -dval(&eps)) Branch (2490:17): [True: 1.84k, False: 7]
|
2491 | goto no_digits; |
2492 | goto fast_failed; |
2493 | } |
2494 | if (leftright) { Branch (2494:13): [True: 0, False: 27.4k]
|
2495 | /* Use Steele & White method of only |
2496 | * generating digits needed. |
2497 | */ |
2498 | dval(&eps) = 0.5/tens[ilim-1] - dval(&eps); |
2499 | for(i = 0;;) { |
2500 | L = (Long)dval(&u); |
2501 | dval(&u) -= L; |
2502 | *s++ = '0' + (int)L; |
2503 | if (dval(&u) < dval(&eps)) Branch (2503:21): [True: 0, False: 0]
|
2504 | goto ret1; |
2505 | if (1. - dval(&u) < dval(&eps)) Branch (2505:21): [True: 0, False: 0]
|
2506 | goto bump_up; |
2507 | if (++i >= ilim) Branch (2507:21): [True: 0, False: 0]
|
2508 | break; |
2509 | dval(&eps) *= 10.; |
2510 | dval(&u) *= 10.; |
2511 | } |
2512 | } |
2513 | else { |
2514 | /* Generate ilim digits, then fix them up. */ |
2515 | dval(&eps) *= tens[ilim-1]; |
2516 | for(i = 1;; i++, 66.1k dval66.1k (&u) *= 10.) { |
2517 | L = (Long)(dval(&u)); |
2518 | if (!(dval(&u) -= L)) Branch (2518:21): [True: 1.83k, False: 91.7k]
|
2519 | ilim = i; |
2520 | *s++ = '0' + (int)L; |
2521 | if (i == ilim) { Branch (2521:21): [True: 27.4k, False: 66.1k]
|
2522 | if (dval(&u) > 0.5 + dval(&eps)) Branch (2522:25): [True: 11.2k, False: 16.1k]
|
2523 | goto bump_up; |
2524 | else if (dval(&u) < 0.5 - dval(&eps)) { Branch (2524:30): [True: 14.9k, False: 1.18k]
|
2525 | while(*--s == '0');2.08k Branch (2525:31): [True: 2.08k, False: 14.9k]
|
2526 | s++; |
2527 | goto ret1; |
2528 | } |
2529 | break; |
2530 | } |
2531 | } |
2532 | } |
2533 | fast_failed: |
2534 | s = s0; |
2535 | dval(&u) = dval(&d2); |
2536 | k = k0; |
2537 | ilim = ilim0; |
2538 | } |
2539 | |
2540 | /* Do we have a "small" integer? */ |
2541 | |
2542 | if (be >= 0 && k <= 67.4k Int_max67.4k ) { Branch (2542:9): [True: 67.4k, False: 1.30M]
Branch (2542:20): [True: 7.33k, False: 60.0k]
|
2543 | /* Yes. */ |
2544 | ds = tens[k]; |
2545 | if (ndigits < 0 && ilim <= 08 ) { Branch (2545:13): [True: 8, False: 7.32k]
Branch (2545:28): [True: 0, False: 8]
|
2546 | S = mhi = 0; |
2547 | if (ilim < 0 || dval(&u) <= 5*ds) Branch (2547:17): [True: 0, False: 0]
Branch (2547:29): [True: 0, False: 0]
|
2548 | goto no_digits; |
2549 | goto one_digit; |
2550 | } |
2551 | for(i = 1;; 7.33k i++, 20.9k dval20.9k (&u) *= 10.) { |
2552 | L = (Long)(dval(&u) / ds); |
2553 | dval(&u) -= L*ds; |
2554 | *s++ = '0' + (int)L; |
2555 | if (!dval(&u)) { Branch (2555:17): [True: 7.31k, False: 20.9k]
|
2556 | break; |
2557 | } |
2558 | if (i == ilim) { Branch (2558:17): [True: 18, False: 20.9k]
|
2559 | dval(&u) += dval(&u); |
2560 | if (dval(&u) > ds || (dval(&u) == ds && L & 1)) { Branch (2560:21): [True: 0, False: 18]
Branch (2560:39): [True: 18, False: 0]
Branch (2560:57): [True: 8, False: 10]
|
2561 | bump_up: |
2562 | while(*--s == '9') Branch (2562:27): [True: 1.21k, False: 11.1k]
|
2563 | if (s == s0) { Branch (2563:29): [True: 96, False: 1.12k]
|
2564 | k++; |
2565 | *s = '0'; |
2566 | break; |
2567 | } |
2568 | ++*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 | while (s > s0 && s[-1] == '0') { Branch (2574:28): [True: 20, False: 0]
Branch (2574:38): [True: 10, False: 10]
|
2575 | --s; |
2576 | } |
2577 | } |
2578 | break; |
2579 | } |
2580 | } |
2581 | goto ret1; |
2582 | } |
2583 | |
2584 | m2 = b2; |
2585 | m5 = b5; |
2586 | if (leftright) { Branch (2586:9): [True: 136k, False: 1.22M]
|
2587 | i = |
2588 | denorm ? be + (989 Bias989 + (P989 -1) - 1 + 1) : Branch (2588:13): [True: 989, False: 135k]
|
2589 | 1 + 135k P135k - bbits; |
2590 | b2 += i; |
2591 | s2 += i; |
2592 | mhi = i2b(1); |
2593 | if (mhi == NULL) Branch (2593:13): [True: 0, False: 136k]
|
2594 | goto failed_malloc; |
2595 | } |
2596 | if (m2 > 0 && s2 > 01.33M ) { Branch (2596:9): [True: 1.33M, False: 28.2k]
Branch (2596:19): [True: 1.33M, False: 0]
|
2597 | i = m2 < s2 ? m21.29M : s244.7k ; Branch (2597:13): [True: 1.29M, False: 44.7k]
|
2598 | b2 -= i; |
2599 | m2 -= i; |
2600 | s2 -= i; |
2601 | } |
2602 | if (b5 > 0) { Branch (2602:9): [True: 1.27M, False: 86.4k]
|
2603 | if (leftright) { Branch (2603:13): [True: 64.5k, False: 1.21M]
|
2604 | if (m5 > 0) { Branch (2604:17): [True: 64.5k, False: 0]
|
2605 | mhi = pow5mult(mhi, m5); |
2606 | if (mhi == NULL) Branch (2606:21): [True: 0, False: 64.5k]
|
2607 | goto failed_malloc; |
2608 | b1 = mult(mhi, b); |
2609 | Bfree(b); |
2610 | b = b1; |
2611 | if (b == NULL) Branch (2611:21): [True: 0, False: 64.5k]
|
2612 | goto failed_malloc; |
2613 | } |
2614 | if ((j = b5 - m5)) { Branch (2614:17): [True: 0, False: 64.5k]
|
2615 | b = pow5mult(b, j); |
2616 | if (b == NULL) Branch (2616:21): [True: 0, False: 0]
|
2617 | goto failed_malloc; |
2618 | } |
2619 | } |
2620 | else { |
2621 | b = pow5mult(b, b5); |
2622 | if (b == NULL) Branch (2622:17): [True: 0, False: 1.21M]
|
2623 | goto failed_malloc; |
2624 | } |
2625 | } |
2626 | S = i2b(1); |
2627 | if (S == NULL) Branch (2627:9): [True: 0, False: 1.36M]
|
2628 | goto failed_malloc; |
2629 | if (s5 > 0) { Branch (2629:9): [True: 80.0k, False: 1.28M]
|
2630 | S = pow5mult(S, s5); |
2631 | if (S == NULL) Branch (2631:13): [True: 0, False: 80.0k]
|
2632 | goto failed_malloc; |
2633 | } |
2634 | |
2635 | /* Check for special case that d is a normalized power of 2. */ |
2636 | |
2637 | spec_case = 0; |
2638 | if ((mode < 2 || leftright1.22M ) Branch (2638:10): [True: 136k, False: 1.22M]
Branch (2638:22): [True: 0, False: 1.22M]
|
2639 | ) { |
2640 | if (!word1(&u) && !(2.53k word02.53k (&u) & Bndry_mask2.53k ) Branch (2640:13): [True: 2.53k, False: 134k]
Branch (2640:27): [True: 2.33k, False: 196]
|
2641 | && word02.33k (&u) & (2.33k Exp_mask2.33k & ~Exp_msk12.33k ) Branch (2641:16): [True: 2.33k, False: 1]
|
2642 | ) { |
2643 | /* The special case */ |
2644 | b2 += Log2P; |
2645 | s2 += Log2P; |
2646 | 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 | i = dshift(S, s2); |
2659 | b2 += i; |
2660 | m2 += i; |
2661 | s2 += i; |
2662 | if (b2 > 0) { Branch (2662:9): [True: 1.36M, False: 26]
|
2663 | b = lshift(b, b2); |
2664 | if (b == NULL) Branch (2664:13): [True: 0, False: 1.36M]
|
2665 | goto failed_malloc; |
2666 | } |
2667 | if (s2 > 0) { Branch (2667:9): [True: 1.36M, False: 1.23k]
|
2668 | S = lshift(S, s2); |
2669 | if (S == NULL) Branch (2669:13): [True: 0, False: 1.36M]
|
2670 | goto failed_malloc; |
2671 | } |
2672 | if (k_check) { Branch (2672:9): [True: 1.32M, False: 42.6k]
|
2673 | if (cmp(b,S) < 0) { Branch (2673:13): [True: 1.60k, False: 1.31M]
|
2674 | k--; |
2675 | b = multadd(b, 10, 0); /* we botched the k estimate */ |
2676 | if (b == NULL) Branch (2676:17): [True: 0, False: 1.60k]
|
2677 | goto failed_malloc; |
2678 | if (leftright) { Branch (2678:17): [True: 1.58k, False: 16]
|
2679 | mhi = multadd(mhi, 10, 0); |
2680 | if (mhi == NULL) Branch (2680:21): [True: 0, False: 1.58k]
|
2681 | goto failed_malloc; |
2682 | } |
2683 | ilim = ilim1; |
2684 | } |
2685 | } |
2686 | if (ilim <= 0 && (1.34M mode == 31.34M || mode == 5136k )) { Branch (2686:9): [True: 1.34M, False: 14.9k]
Branch (2686:23): [True: 1.21M, False: 136k]
Branch (2686:36): [True: 0, False: 136k]
|
2687 | if (ilim < 0) { Branch (2687:13): [True: 1.21M, False: 8]
|
2688 | /* no digits, fcvt style */ |
2689 | no_digits: |
2690 | k = -1 - ndigits; |
2691 | goto ret; |
2692 | } |
2693 | else { |
2694 | S = multadd(S, 5, 0); |
2695 | if (S == NULL) Branch (2695:17): [True: 0, False: 8]
|
2696 | goto failed_malloc; |
2697 | if (cmp(b, S) <= 0) Branch (2697:17): [True: 1, False: 7]
|
2698 | goto no_digits; |
2699 | } |
2700 | one_digit: |
2701 | *s++ = '1'; |
2702 | k++; |
2703 | goto ret; |
2704 | } |
2705 | if (leftright) { Branch (2705:9): [True: 136k, False: 14.9k]
|
2706 | if (m2 > 0) { Branch (2706:13): [True: 135k, False: 1.57k]
|
2707 | mhi = lshift(mhi, m2); |
2708 | if (mhi == NULL) Branch (2708:17): [True: 0, False: 135k]
|
2709 | goto failed_malloc; |
2710 | } |
2711 | |
2712 | /* Compute mlo -- check for special case |
2713 | * that d is a normalized power of 2. |
2714 | */ |
2715 | |
2716 | mlo = mhi; |
2717 | if (spec_case) { Branch (2717:13): [True: 2.33k, False: 134k]
|
2718 | mhi = Balloc(mhi->k); |
2719 | if (mhi == NULL) Branch (2719:17): [True: 0, False: 2.33k]
|
2720 | goto failed_malloc; |
2721 | Bcopy(mhi, mlo); |
2722 | mhi = lshift(mhi, Log2P); |
2723 | if (mhi == NULL) Branch (2723:17): [True: 0, False: 2.33k]
|
2724 | goto failed_malloc; |
2725 | } |
2726 | |
2727 | for(i = 1;;136k i++1.96M ) { |
2728 | dig = quorem(b,S) + '0'; |
2729 | /* Do we yet have the shortest decimal string |
2730 | * that will round to d? |
2731 | */ |
2732 | j = cmp(b, mlo); |
2733 | delta = diff(S, mhi); |
2734 | if (delta == NULL) Branch (2734:17): [True: 0, False: 2.09M]
|
2735 | goto failed_malloc; |
2736 | j1 = delta->sign ? 135.2k : cmp(b, delta)2.06M ; Branch (2736:18): [True: 35.2k, False: 2.06M]
|
2737 | Bfree(delta); |
2738 | if (j1 == 0 && mode != 12.14k && !(2.14k word12.14k (&u) & 1) Branch (2738:17): [True: 2.14k, False: 2.09M]
Branch (2738:28): [True: 2.14k, False: 0]
Branch (2738:41): [True: 1.41k, False: 722]
|
2739 | ) { |
2740 | if (dig == '9') Branch (2740:21): [True: 1, False: 1.41k]
|
2741 | goto round_9_up; |
2742 | if (j > 0) Branch (2742:21): [True: 236, False: 1.18k]
|
2743 | dig++; |
2744 | *s++ = dig; |
2745 | goto ret; |
2746 | } |
2747 | if (j < 0 || (2.00M j == 02.00M && mode != 1739 Branch (2747:17): [True: 95.1k, False: 2.00M]
Branch (2747:27): [True: 739, False: 2.00M]
Branch (2747:37): [True: 739, False: 0]
|
2748 | && !(739 word1739 (&u) & 1) Branch (2748:30): [True: 511, False: 228]
|
2749 | )) { |
2750 | if (!b->x[0] && b->wds <= 125.0k ) { Branch (2750:21): [True: 25.0k, False: 70.5k]
Branch (2750:33): [True: 6.78k, False: 18.3k]
|
2751 | goto accept_dig; |
2752 | } |
2753 | if (j1 > 0) { Branch (2753:21): [True: 45.6k, False: 43.1k]
|
2754 | b = lshift(b, 1); |
2755 | if (b == NULL) Branch (2755:25): [True: 0, False: 45.6k]
|
2756 | goto failed_malloc; |
2757 | j1 = cmp(b, S); |
2758 | if ((j1 > 0 || (23.3k j1 == 023.3k && dig & 170 )) Branch (2758:26): [True: 22.3k, False: 23.3k]
Branch (2758:37): [True: 70, False: 23.2k]
Branch (2758:48): [True: 34, False: 36]
|
2759 | && dig++ == '9'22.3k ) Branch (2759:28): [True: 96, False: 22.2k]
|
2760 | goto round_9_up; |
2761 | } |
2762 | accept_dig: |
2763 | *s++ = dig; |
2764 | goto ret; |
2765 | } |
2766 | if (j1 > 0) { Branch (2766:17): [True: 39.9k, False: 1.96M]
|
2767 | if (dig == '9') { /* possible if i == 1 */ Branch (2767:21): [True: 666, False: 39.2k]
|
2768 | round_9_up: |
2769 | *s++ = '9'; |
2770 | goto roundoff; |
2771 | } |
2772 | *s++ = dig + 1; |
2773 | goto ret; |
2774 | } |
2775 | *s++ = dig; |
2776 | if (i == ilim) Branch (2776:17): [True: 0, False: 1.96M]
|
2777 | break; |
2778 | b = multadd(b, 10, 0); |
2779 | if (b == NULL) Branch (2779:17): [True: 0, False: 1.96M]
|
2780 | goto failed_malloc; |
2781 | if (mlo == mhi) { Branch (2781:17): [True: 1.93M, False: 32.1k]
|
2782 | mlo = mhi = multadd(mhi, 10, 0); |
2783 | if (mlo == NULL) Branch (2783:21): [True: 0, False: 1.93M]
|
2784 | goto failed_malloc; |
2785 | } |
2786 | else { |
2787 | mlo = multadd(mlo, 10, 0); |
2788 | if (mlo == NULL) Branch (2788:21): [True: 0, False: 32.1k]
|
2789 | goto failed_malloc; |
2790 | mhi = multadd(mhi, 10, 0); |
2791 | if (mhi == NULL) Branch (2791:21): [True: 0, False: 32.1k]
|
2792 | goto failed_malloc; |
2793 | } |
2794 | } |
2795 | } |
2796 | else |
2797 | for(i = 1;; 14.9k i++489k ) { |
2798 | *s++ = dig = quorem(b,S) + '0'; |
2799 | if (!b->x[0] && b->wds <= 182.0k ) { Branch (2799:17): [True: 82.0k, False: 422k]
Branch (2799:29): [True: 11.6k, False: 70.4k]
|
2800 | goto ret; |
2801 | } |
2802 | if (i >= ilim) Branch (2802:17): [True: 3.31k, False: 489k]
|
2803 | break; |
2804 | b = multadd(b, 10, 0); |
2805 | if (b == NULL) Branch (2805:17): [True: 0, False: 489k]
|
2806 | goto failed_malloc; |
2807 | } |
2808 | |
2809 | /* Round off last digit */ |
2810 | |
2811 | b = lshift(b, 1); |
2812 | if (b == NULL) Branch (2812:9): [True: 0, False: 3.31k]
|
2813 | goto failed_malloc; |
2814 | j = cmp(b, S); |
2815 | if (j > 0 || (1.24k j == 01.24k && dig & 1140 )) { Branch (2815:9): [True: 2.06k, False: 1.24k]
Branch (2815:19): [True: 140, False: 1.10k]
Branch (2815:29): [True: 78, False: 62]
|
2816 | roundoff: |
2817 | while(*--s == '9') Branch (2817:15): [True: 913, False: 2.13k]
|
2818 | if (s == s0) { Branch (2818:17): [True: 771, False: 142]
|
2819 | k++; |
2820 | *s++ = '1'; |
2821 | goto ret; |
2822 | } |
2823 | ++*s++; |
2824 | } |
2825 | else { |
2826 | while(*--s == '0');188 Branch (2826:15): [True: 188, False: 1.17k]
|
2827 | s++; |
2828 | } |
2829 | ret: |
2830 | Bfree(S); |
2831 | if (mhi) { Branch (2831:9): [True: 136k, False: 1.22M]
|
2832 | if (mlo && mlo != mhi) Branch (2832:13): [True: 136k, False: 0]
Branch (2832:20): [True: 2.33k, False: 134k]
|
2833 | Bfree(mlo); |
2834 | Bfree(mhi); |
2835 | } |
2836 | ret1: |
2837 | Bfree(b); |
2838 | *s = 0; |
2839 | *decpt = k + 1; |
2840 | if (rve) Branch (2840:9): [True: 1.39M, False: 0]
|
2841 | *rve = s; |
2842 | return s0; |
2843 | failed_malloc: |
2844 | if (S) Branch (2844:9): [True: 0, False: 0]
|
2845 | Bfree(S); |
2846 | if (mlo && mlo != mhi) Branch (2846:9): [True: 0, False: 0]
Branch (2846:16): [True: 0, False: 0]
|
2847 | Bfree(mlo); |
2848 | if (mhi) Branch (2848:9): [True: 0, False: 0]
|
2849 | Bfree(mhi); |
2850 | if (b) Branch (2850:9): [True: 0, False: 0]
|
2851 | Bfree(b); |
2852 | if (s0) Branch (2852:9): [True: 0, False: 0]
|
2853 | _Py_dg_freedtoa(s0); |
2854 | return NULL; |
2855 | } |
2856 | #ifdef __cplusplus |
2857 | } |
2858 | #endif |
2859 | |
2860 | #endif // _PY_SHORT_FLOAT_REPR == 1 |