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