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