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