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