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