Python code coverage for Python/dtoa.c

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