1	/****************************************************************
2	*
3	* The author of this software is David M. Gay.
4	*
5	* Copyright (c) 1991, 2000, 2001 by Lucent Technologies.
6	* Copyright (C) 2002, 2005, 2006, 2007, 2008 Apple Inc. All rights reserved.
7	*
8	* Permission to use, copy, modify, and distribute this software for any
9	* purpose without fee is hereby granted, provided that this entire notice
10	* is included in all copies of any software which is or includes a copy
11	* or modification of this software and in all copies of the supporting
12	* documentation for such software.
13	*
14	* THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
15	* WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY
16	* REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
17	* OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
18	*
19	***************************************************************/
20
21	/* Please send bug reports to
22	David M. Gay
23	Bell Laboratories, Room 2C-463
24	600 Mountain Avenue
25	Murray Hill, NJ 07974-0636
26	U.S.A.
27	dmg@bell-labs.com
28	*/
29
30	/* On a machine with IEEE extended-precision registers, it is
31	* necessary to specify double-precision (53-bit) rounding precision
32	* before invoking strtod or dtoa. If the machine uses (the equivalent
33	* of) Intel 80x87 arithmetic, the call
34	* _control87(PC_53, MCW_PC);
35	* does this with many compilers. Whether this or another call is
36	* appropriate depends on the compiler; for this to work, it may be
37	* necessary to #include "float.h" or another system-dependent header
38	* file.
39	*/
40
41	/* strtod for IEEE-arithmetic machines.
42	*
43	* This strtod returns a nearest machine number to the input decimal
44	* string (or sets errno to ERANGE). With IEEE arithmetic, ties are
45	* broken by the IEEE round-even rule. Otherwise ties are broken by
46	* biased rounding (add half and chop).
47	*
48	* Inspired loosely by William D. Clinger's paper "How to Read Floating
49	* Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
50	*
51	* Modifications:
52	*
53	* 1. We only require IEEE.
54	* 2. We get by with floating-point arithmetic in a case that
55	* Clinger missed -- when we're computing d * 10^n
56	* for a small integer d and the integer n is not too
57	* much larger than 22 (the maximum integer k for which
58	* we can represent 10^k exactly), we may be able to
59	* compute (d*10^k) * 10^(e-k) with just one roundoff.
60	* 3. Rather than a bit-at-a-time adjustment of the binary
61	* result in the hard case, we use floating-point
62	* arithmetic to determine the adjustment to within
63	* one bit; only in really hard cases do we need to
64	* compute a second residual.
65	* 4. Because of 3., we don't need a large table of powers of 10
66	* for ten-to-e (just some small tables, e.g. of 10^k
67	* for 0 <= k <= 22).
68	*/
69
70	/*
71	* #define IEEE_8087 for IEEE-arithmetic machines where the least
72	* significant byte has the lowest address.
73	* #define IEEE_MC68k for IEEE-arithmetic machines where the most
74	* significant byte has the lowest address.
75	* #define No_leftright to omit left-right logic in fast floating-point
76	* computation of dtoa.
77	* #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
78	* and Honor_FLT_ROUNDS is not #defined.
79	* #define Inaccurate_Divide for IEEE-format with correctly rounded
80	* products but inaccurate quotients, e.g., for Intel i860.
81	* #define USE_LONG_LONG on machines that have a "long long"
82	* integer type (of >= 64 bits), and performance testing shows that
83	* it is faster than 32-bit fallback (which is often not the case
84	* on 32-bit machines). On such machines, you can #define Just_16
85	* to store 16 bits per 32-bit int32_t when doing high-precision integer
86	* arithmetic. Whether this speeds things up or slows things down
87	* depends on the machine and the number being converted.
88	* #define Bad_float_h if your system lacks a float.h or if it does not
89	* define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
90	* FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
91	* #define INFNAN_CHECK on IEEE systems to cause strtod to check for
92	* Infinity and NaN (case insensitively). On some systems (e.g.,
93	* some HP systems), it may be necessary to #define NAN_WORD0
94	* appropriately -- to the most significant word of a quiet NaN.
95	* (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.)
96	* When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined,
97	* strtod also accepts (case insensitively) strings of the form
98	* NaN(x), where x is a string of hexadecimal digits and spaces;
99	* if there is only one string of hexadecimal digits, it is taken
100	* for the 52 fraction bits of the resulting NaN; if there are two
101	* or more strings of hex digits, the first is for the high 20 bits,
102	* the second and subsequent for the low 32 bits, with intervening
103	* white space ignored; but if this results in none of the 52
104	* fraction bits being on (an IEEE Infinity symbol), then NAN_WORD0
105	* and NAN_WORD1 are used instead.
106	* #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that
107	* avoids underflows on inputs whose result does not underflow.
108	* If you #define NO_IEEE_Scale on a machine that uses IEEE-format
109	* floating-point numbers and flushes underflows to zero rather
110	* than implementing gradual underflow, then you must also #define
111	* Sudden_Underflow.
112	* #define YES_ALIAS to permit aliasing certain double values with
113	* arrays of ULongs. This leads to slightly better code with
114	* some compilers and was always used prior to 19990916, but it
115	* is not strictly legal and can cause trouble with aggressively
116	* optimizing compilers (e.g., gcc 2.95.1 under -O2).
117	* #define SET_INEXACT if IEEE arithmetic is being used and extra
118	* computation should be done to set the inexact flag when the
119	* result is inexact and avoid setting inexact when the result
120	* is exact. In this case, dtoa.c must be compiled in
121	* an environment, perhaps provided by #include "dtoa.c" in a
122	* suitable wrapper, that defines two functions,
123	* int get_inexact(void);
124	* void clear_inexact(void);
125	* such that get_inexact() returns a nonzero value if the
126	* inexact bit is already set, and clear_inexact() sets the
127	* inexact bit to 0. When SET_INEXACT is #defined, strtod
128	* also does extra computations to set the underflow and overflow
129	* flags when appropriate (i.e., when the result is tiny and
130	* inexact or when it is a numeric value rounded to +-infinity).
131	* #define NO_ERRNO if strtod should not assign errno = ERANGE when
132	* the result overflows to +-Infinity or underflows to 0.
133	*/
134
135	#include "config.h"
136	#include "dtoa.h"
137
138	#if HAVE(ERRNO_H)
139	#include <errno.h>
140	#else
141	#define NO_ERRNO
142	#endif
143	#include <math.h>
144	#include <stdint.h>
145	#include <stdlib.h>
146	#include <string.h>
147	#include <wtf/AlwaysInline.h>
148	#include <wtf/Assertions.h>
149	#include <wtf/FastMalloc.h>
150	#include <wtf/Vector.h>
151	#include <wtf/Threading.h>
152
153	#include <stdio.h>
154
155	#include <wtf/MathExtras.h>
156
157	#if COMPILER(MSVC)
158	#pragma warning(disable: 4244)
159	#pragma warning(disable: 4245)
160	#pragma warning(disable: 4554)
161	#endif
162
163	#if CPU(BIG_ENDIAN)
164	#define IEEE_MC68k
165	#elif CPU(MIDDLE_ENDIAN)
166	#define IEEE_ARM
167	#else
168	#define IEEE_8087
169	#endif
170
171	#define INFNAN_CHECK
172
173	#if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(IEEE_ARM) != 1
174	Exactly one of IEEE_8087, IEEE_ARM or IEEE_MC68k should be defined.
175	#endif
176
177	namespace WTF {
178
179	#if ENABLE(JSC_MULTIPLE_THREADS)
180	Mutex* s_dtoaP5Mutex;
181	#endif
182
183	typedef union { double d; uint32_t L[2]; } U;
184
185	#ifdef YES_ALIAS
186	#define dval(x) x
187	#ifdef IEEE_8087
188	#define word0(x) ((uint32_t*)&x)[1]
189	#define word1(x) ((uint32_t*)&x)[0]
190	#else
191	#define word0(x) ((uint32_t*)&x)[0]
192	#define word1(x) ((uint32_t*)&x)[1]
193	#endif
194	#else
195	#ifdef IEEE_8087
196	#define word0(x) (x)->L[1]
197	#define word1(x) (x)->L[0]
198	#else
199	#define word0(x) (x)->L[0]
200	#define word1(x) (x)->L[1]
201	#endif
202	#define dval(x) (x)->d
203	#endif
204
205	/* The following definition of Storeinc is appropriate for MIPS processors.
206	* An alternative that might be better on some machines is
207	* #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
208	*/
209	#if defined(IEEE_8087) || defined(IEEE_ARM)
210	#define Storeinc(a,b,c) (((unsigned short*)a)[1] = (unsigned short)b, ((unsigned short*)a)[0] = (unsigned short)c, a++)
211	#else
212	#define Storeinc(a,b,c) (((unsigned short*)a)[0] = (unsigned short)b, ((unsigned short*)a)[1] = (unsigned short)c, a++)
213	#endif
214
215	#define Exp_shift 20
216	#define Exp_shift1 20
217	#define Exp_msk1 0x100000
218	#define Exp_msk11 0x100000
219	#define Exp_mask 0x7ff00000
220	#define P 53
221	#define Bias 1023
222	#define Emin (-1022)
223	#define Exp_1 0x3ff00000
224	#define Exp_11 0x3ff00000
225	#define Ebits 11
226	#define Frac_mask 0xfffff
227	#define Frac_mask1 0xfffff
228	#define Ten_pmax 22
229	#define Bletch 0x10
230	#define Bndry_mask 0xfffff
231	#define Bndry_mask1 0xfffff
232	#define LSB 1
233	#define Sign_bit 0x80000000
234	#define Log2P 1
235	#define Tiny0 0
236	#define Tiny1 1
237	#define Quick_max 14
238	#define Int_max 14
239
240	#if !defined(NO_IEEE_Scale)
241	#undef Avoid_Underflow
242	#define Avoid_Underflow
243	#endif
244
245	#if !defined(Flt_Rounds)
246	#if defined(FLT_ROUNDS)
247	#define Flt_Rounds FLT_ROUNDS
248	#else
249	#define Flt_Rounds 1
250	#endif
251	#endif /*Flt_Rounds*/
252
253
254	#define rounded_product(a,b) a *= b
255	#define rounded_quotient(a,b) a /= b
256
257	#define Big0 (Frac_mask1 | Exp_msk1 * (DBL_MAX_EXP + Bias - 1))
258	#define Big1 0xffffffff
259
260
261	// FIXME: we should remove non-Pack_32 mode since it is unused and unmaintained
262	#ifndef Pack_32
263	#define Pack_32
264	#endif
265
266	#if CPU(PPC64) || CPU(X86_64)
267	// FIXME: should we enable this on all 64-bit CPUs?
268	// 64-bit emulation provided by the compiler is likely to be slower than dtoa own code on 32-bit hardware.
269	#define USE_LONG_LONG
270	#endif
271
272	#ifndef USE_LONG_LONG
273	#ifdef Just_16
274	#undef Pack_32
275	/* When Pack_32 is not defined, we store 16 bits per 32-bit int32_t.
276	* This makes some inner loops simpler and sometimes saves work
277	* during multiplications, but it often seems to make things slightly
278	* slower. Hence the default is now to store 32 bits per int32_t.
279	*/
280	#endif
281	#endif
282
283	#define Kmax 15
284
285	struct BigInt {
286	BigInt() : sign(0) { }
287	int sign;
288
289	void clear()
290	{
291	sign = 0;
292	m_words.clear();
293	}
294
295	size_t size() const
296	{
297	return m_words.size();
298	}
299
300	void resize(size_t s)
301	{
302	m_words.resize(size: s);
303	}
304
305	uint32_t* words()
306	{
307	return m_words.data();
308	}
309
310	const uint32_t* words() const
311	{
312	return m_words.data();
313	}
314
315	void append(uint32_t w)
316	{
317	m_words.append(val: w);
318	}
319
320	Vector<uint32_t, 16> m_words;
321	};
322
323	static void multadd(BigInt& b, int m, int a) /* multiply by m and add a */
324	{
325	#ifdef USE_LONG_LONG
326	unsigned long long carry;
327	#else
328	uint32_t carry;
329	#endif
330
331	int wds = b.size();
332	uint32_t* x = b.words();
333	int i = 0;
334	carry = a;
335	do {
336	#ifdef USE_LONG_LONG
337	unsigned long long y = *x * (unsigned long long)m + carry;
338	carry = y >> 32;
339	*x++ = (uint32_t)y & 0xffffffffUL;
340	#else
341	#ifdef Pack_32
342	uint32_t xi = *x;
343	uint32_t y = (xi & 0xffff) * m + carry;
344	uint32_t z = (xi >> 16) * m + (y >> 16);
345	carry = z >> 16;
346	*x++ = (z << 16) + (y & 0xffff);
347	#else
348	uint32_t y = *x * m + carry;
349	carry = y >> 16;
350	*x++ = y & 0xffff;
351	#endif
352	#endif
353	} while (++i < wds);
354
355	if (carry)
356	b.append(w: (uint32_t)carry);
357	}
358
359	static void s2b(BigInt& b, const char* s, int nd0, int nd, uint32_t y9)
360	{
361	int k;
362	int32_t y;
363	int32_t x = (nd + 8) / 9;
364
365	for (k = 0, y = 1; x > y; y <<= 1, k++) { }
366	#ifdef Pack_32
367	b.sign = 0;
368	b.resize(s: 1);
369	b.words()[0] = y9;
370	#else
371	b.sign = 0;
372	b.resize((b->x[1] = y9 >> 16) ? 2 : 1);
373	b.words()[0] = y9 & 0xffff;
374	#endif
375
376	int i = 9;
377	if (9 < nd0) {
378	s += 9;
379	do {
380	multadd(b, m: 10, a: *s++ - '0');
381	} while (++i < nd0);
382	s++;
383	} else
384	s += 10;
385	for (; i < nd; i++)
386	multadd(b, m: 10, a: *s++ - '0');
387	}
388
389	static int hi0bits(uint32_t x)
390	{
391	int k = 0;
392
393	if (!(x & 0xffff0000)) {
394	k = 16;
395	x <<= 16;
396	}
397	if (!(x & 0xff000000)) {
398	k += 8;
399	x <<= 8;
400	}
401	if (!(x & 0xf0000000)) {
402	k += 4;
403	x <<= 4;
404	}
405	if (!(x & 0xc0000000)) {
406	k += 2;
407	x <<= 2;
408	}
409	if (!(x & 0x80000000)) {
410	k++;
411	if (!(x & 0x40000000))
412	return 32;
413	}
414	return k;
415	}
416
417	static int lo0bits (uint32_t* y)
418	{
419	int k;
420	uint32_t x = *y;
421
422	if (x & 7) {
423	if (x & 1)
424	return 0;
425	if (x & 2) {
426	*y = x >> 1;
427	return 1;
428	}
429	*y = x >> 2;
430	return 2;
431	}
432	k = 0;
433	if (!(x & 0xffff)) {
434	k = 16;
435	x >>= 16;
436	}
437	if (!(x & 0xff)) {
438	k += 8;
439	x >>= 8;
440	}
441	if (!(x & 0xf)) {
442	k += 4;
443	x >>= 4;
444	}
445	if (!(x & 0x3)) {
446	k += 2;
447	x >>= 2;
448	}
449	if (!(x & 1)) {
450	k++;
451	x >>= 1;
452	if (!x & 1)
453	return 32;
454	}
455	*y = x;
456	return k;
457	}
458
459	static void i2b(BigInt& b, int i)
460	{
461	b.sign = 0;
462	b.resize(s: 1);
463	b.words()[0] = i;
464	}
465
466	static void mult(BigInt& aRef, const BigInt& bRef)
467	{
468	const BigInt* a = &aRef;
469	const BigInt* b = &bRef;
470	BigInt c;
471	int wa, wb, wc;
472	const uint32_t *x = 0, *xa, *xb, *xae, *xbe;
473	uint32_t *xc, *xc0;
474	uint32_t y;
475	#ifdef USE_LONG_LONG
476	unsigned long long carry, z;
477	#else
478	uint32_t carry, z;
479	#endif
480
481	if (a->size() < b->size()) {
482	const BigInt* tmp = a;
483	a = b;
484	b = tmp;
485	}
486
487	wa = a->size();
488	wb = b->size();
489	wc = wa + wb;
490	c.resize(s: wc);
491
492	for (xc = c.words(), xa = xc + wc; xc < xa; xc++)
493	*xc = 0;
494	xa = a->words();
495	xae = xa + wa;
496	xb = b->words();
497	xbe = xb + wb;
498	xc0 = c.words();
499	#ifdef USE_LONG_LONG
500	for (; xb < xbe; xc0++) {
501	if ((y = *xb++)) {
502	x = xa;
503	xc = xc0;
504	carry = 0;
505	do {
506	z = *x++ * (unsigned long long)y + *xc + carry;
507	carry = z >> 32;
508	*xc++ = (uint32_t)z & 0xffffffffUL;
509	} while (x < xae);
510	*xc = (uint32_t)carry;
511	}
512	}
513	#else
514	#ifdef Pack_32
515	for (; xb < xbe; xb++, xc0++) {
516	if ((y = *xb & 0xffff)) {
517	x = xa;
518	xc = xc0;
519	carry = 0;
520	do {
521	z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
522	carry = z >> 16;
523	uint32_t z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
524	carry = z2 >> 16;
525	Storeinc(xc, z2, z);
526	} while (x < xae);
527	*xc = carry;
528	}
529	if ((y = *xb >> 16)) {
530	x = xa;
531	xc = xc0;
532	carry = 0;
533	uint32_t z2 = *xc;
534	do {
535	z = (*x & 0xffff) * y + (*xc >> 16) + carry;
536	carry = z >> 16;
537	Storeinc(xc, z, z2);
538	z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
539	carry = z2 >> 16;
540	} while (x < xae);
541	*xc = z2;
542	}
543	}
544	#else
545	for(; xb < xbe; xc0++) {
546	if ((y = *xb++)) {
547	x = xa;
548	xc = xc0;
549	carry = 0;
550	do {
551	z = *x++ * y + *xc + carry;
552	carry = z >> 16;
553	*xc++ = z & 0xffff;
554	} while (x < xae);
555	*xc = carry;
556	}
557	}
558	#endif
559	#endif
560	for (xc0 = c.words(), xc = xc0 + wc; wc > 0 && !*--xc; --wc) { }
561	c.resize(s: wc);
562	aRef = c;
563	}
564
565	struct P5Node : Noncopyable {
566	BigInt val;
567	P5Node* next;
568	};
569
570	static P5Node* p5s;
571	static int p5s_count;
572
573	static ALWAYS_INLINE void pow5mult(BigInt& b, int k)
574	{
575	static int p05[3] = { 5, 25, 125 };
576
577	if (int i = k & 3)
578	multadd(b, m: p05[i - 1], a: 0);
579
580	if (!(k >>= 2))
581	return;
582
583	#if ENABLE(JSC_MULTIPLE_THREADS)
584	s_dtoaP5Mutex->lock();
585	#endif
586	P5Node* p5 = p5s;
587
588	if (!p5) {
589	/* first time */
590	p5 = new P5Node;
591	i2b(b&: p5->val, i: 625);
592	p5->next = 0;
593	p5s = p5;
594	p5s_count = 1;
595	}
596
597	int p5s_count_local = p5s_count;
598	#if ENABLE(JSC_MULTIPLE_THREADS)
599	s_dtoaP5Mutex->unlock();
600	#endif
601	int p5s_used = 0;
602
603	for (;;) {
604	if (k & 1)
605	mult(aRef&: b, bRef: p5->val);
606
607	if (!(k >>= 1))
608	break;
609
610	if (++p5s_used == p5s_count_local) {
611	#if ENABLE(JSC_MULTIPLE_THREADS)
612	s_dtoaP5Mutex->lock();
613	#endif
614	if (p5s_used == p5s_count) {
615	ASSERT(!p5->next);
616	p5->next = new P5Node;
617	p5->next->next = 0;
618	p5->next->val = p5->val;
619	mult(aRef&: p5->next->val, bRef: p5->next->val);
620	++p5s_count;
621	}
622
623	p5s_count_local = p5s_count;
624	#if ENABLE(JSC_MULTIPLE_THREADS)
625	s_dtoaP5Mutex->unlock();
626	#endif
627	}
628	p5 = p5->next;
629	}
630	}
631
632	static ALWAYS_INLINE void lshift(BigInt& b, int k)
633	{
634	#ifdef Pack_32
635	int n = k >> 5;
636	#else
637	int n = k >> 4;
638	#endif
639
640	int origSize = b.size();
641	int n1 = n + origSize + 1;
642
643	if (k &= 0x1f)
644	b.resize(s: b.size() + n + 1);
645	else
646	b.resize(s: b.size() + n);
647
648	const uint32_t* srcStart = b.words();
649	uint32_t* dstStart = b.words();
650	const uint32_t* src = srcStart + origSize - 1;
651	uint32_t* dst = dstStart + n1 - 1;
652	#ifdef Pack_32
653	if (k) {
654	uint32_t hiSubword = 0;
655	int s = 32 - k;
656	for (; src >= srcStart; --src) {
657	*dst-- = hiSubword | *src >> s;
658	hiSubword = *src << k;
659	}
660	*dst = hiSubword;
661	ASSERT(dst == dstStart + n);
662
663	b.resize(s: origSize + n + (b.words()[n1 - 1] != 0));
664	}
665	#else
666	if (k &= 0xf) {
667	uint32_t hiSubword = 0;
668	int s = 16 - k;
669	for (; src >= srcStart; --src) {
670	*dst-- = hiSubword | *src >> s;
671	hiSubword = (*src << k) & 0xffff;
672	}
673	*dst = hiSubword;
674	ASSERT(dst == dstStart + n);
675	result->wds = b->wds + n + (result->x[n1 - 1] != 0);
676	}
677	#endif
678	else {
679	do {
680	*--dst = *src--;
681	} while (src >= srcStart);
682	}
683	for (dst = dstStart + n; dst != dstStart; )
684	*--dst = 0;
685
686	ASSERT(b.size() <= 1 || b.words()[b.size() - 1]);
687	}
688
689	static int cmp(const BigInt& a, const BigInt& b)
690	{
691	const uint32_t *xa, *xa0, *xb, *xb0;
692	int i, j;
693
694	i = a.size();
695	j = b.size();
696	ASSERT(i <= 1 || a.words()[i - 1]);
697	ASSERT(j <= 1 || b.words()[j - 1]);
698	if (i -= j)
699	return i;
700	xa0 = a.words();
701	xa = xa0 + j;
702	xb0 = b.words();
703	xb = xb0 + j;
704	for (;;) {
705	if (*--xa != *--xb)
706	return *xa < *xb ? -1 : 1;
707	if (xa <= xa0)
708	break;
709	}
710	return 0;
711	}
712
713	static ALWAYS_INLINE void diff(BigInt& c, const BigInt& aRef, const BigInt& bRef)
714	{
715	const BigInt* a = &aRef;
716	const BigInt* b = &bRef;
717	int i, wa, wb;
718	uint32_t *xc;
719
720	i = cmp(a: *a, b: *b);
721	if (!i) {
722	c.sign = 0;
723	c.resize(s: 1);
724	c.words()[0] = 0;
725	return;
726	}
727	if (i < 0) {
728	const BigInt* tmp = a;
729	a = b;
730	b = tmp;
731	i = 1;
732	} else
733	i = 0;
734
735	wa = a->size();
736	const uint32_t* xa = a->words();
737	const uint32_t* xae = xa + wa;
738	wb = b->size();
739	const uint32_t* xb = b->words();
740	const uint32_t* xbe = xb + wb;
741
742	c.resize(s: wa);
743	c.sign = i;
744	xc = c.words();
745	#ifdef USE_LONG_LONG
746	unsigned long long borrow = 0;
747	do {
748	unsigned long long y = (unsigned long long)*xa++ - *xb++ - borrow;
749	borrow = y >> 32 & (uint32_t)1;
750	*xc++ = (uint32_t)y & 0xffffffffUL;
751	} while (xb < xbe);
752	while (xa < xae) {
753	unsigned long long y = *xa++ - borrow;
754	borrow = y >> 32 & (uint32_t)1;
755	*xc++ = (uint32_t)y & 0xffffffffUL;
756	}
757	#else
758	uint32_t borrow = 0;
759	#ifdef Pack_32
760	do {
761	uint32_t y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
762	borrow = (y & 0x10000) >> 16;
763	uint32_t z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
764	borrow = (z & 0x10000) >> 16;
765	Storeinc(xc, z, y);
766	} while (xb < xbe);
767	while (xa < xae) {
768	uint32_t y = (*xa & 0xffff) - borrow;
769	borrow = (y & 0x10000) >> 16;
770	uint32_t z = (*xa++ >> 16) - borrow;
771	borrow = (z & 0x10000) >> 16;
772	Storeinc(xc, z, y);
773	}
774	#else
775	do {
776	uint32_t y = *xa++ - *xb++ - borrow;
777	borrow = (y & 0x10000) >> 16;
778	*xc++ = y & 0xffff;
779	} while (xb < xbe);
780	while (xa < xae) {
781	uint32_t y = *xa++ - borrow;
782	borrow = (y & 0x10000) >> 16;
783	*xc++ = y & 0xffff;
784	}
785	#endif
786	#endif
787	while (!*--xc)
788	wa--;
789	c.resize(s: wa);
790	}
791
792	static double ulp(U *x)
793	{
794	int32_t L;
795	U u;
796
797	L = (word0(x) & Exp_mask) - (P - 1) * Exp_msk1;
798	#ifndef Avoid_Underflow
799	#ifndef Sudden_Underflow
800	if (L > 0) {
801	#endif
802	#endif
803	word0(&u) = L;
804	word1(&u) = 0;
805	#ifndef Avoid_Underflow
806	#ifndef Sudden_Underflow
807	} else {
808	L = -L >> Exp_shift;
809	if (L < Exp_shift) {
810	word0(&u) = 0x80000 >> L;
811	word1(&u) = 0;
812	} else {
813	word0(&u) = 0;
814	L -= Exp_shift;
815	word1(&u) = L >= 31 ? 1 : 1 << 31 - L;
816	}
817	}
818	#endif
819	#endif
820	return dval(&u);
821	}
822
823	static double b2d(const BigInt& a, int* e)
824	{
825	const uint32_t* xa;
826	const uint32_t* xa0;
827	uint32_t w;
828	uint32_t y;
829	uint32_t z;
830	int k;
831	U d;
832
833	#define d0 word0(&d)
834	#define d1 word1(&d)
835
836	xa0 = a.words();
837	xa = xa0 + a.size();
838	y = *--xa;
839	ASSERT(y);
840	k = hi0bits(x: y);
841	*e = 32 - k;
842	#ifdef Pack_32
843	if (k < Ebits) {
844	d0 = Exp_1 | (y >> (Ebits - k));
845	w = xa > xa0 ? *--xa : 0;
846	d1 = (y << (32 - Ebits + k)) | (w >> (Ebits - k));
847	goto ret_d;
848	}
849	z = xa > xa0 ? *--xa : 0;
850	if (k -= Ebits) {
851	d0 = Exp_1 | (y << k) | (z >> (32 - k));
852	y = xa > xa0 ? *--xa : 0;
853	d1 = (z << k) | (y >> (32 - k));
854	} else {
855	d0 = Exp_1 | y;
856	d1 = z;
857	}
858	#else
859	if (k < Ebits + 16) {
860	z = xa > xa0 ? *--xa : 0;
861	d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
862	w = xa > xa0 ? *--xa : 0;
863	y = xa > xa0 ? *--xa : 0;
864	d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
865	goto ret_d;
866	}
867	z = xa > xa0 ? *--xa : 0;
868	w = xa > xa0 ? *--xa : 0;
869	k -= Ebits + 16;
870	d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
871	y = xa > xa0 ? *--xa : 0;
872	d1 = w << k + 16 | y << k;
873	#endif
874	ret_d:
875	#undef d0
876	#undef d1
877	return dval(&d);
878	}
879
880	static ALWAYS_INLINE void d2b(BigInt& b, U* d, int* e, int* bits)
881	{
882	int de, k;
883	uint32_t *x, y, z;
884	#ifndef Sudden_Underflow
885	int i;
886	#endif
887	#define d0 word0(d)
888	#define d1 word1(d)
889
890	b.sign = 0;
891	#ifdef Pack_32
892	b.resize(s: 1);
893	#else
894	b.resize(2);
895	#endif
896	x = b.words();
897
898	z = d0 & Frac_mask;
899	d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
900	#ifdef Sudden_Underflow
901	de = (int)(d0 >> Exp_shift);
902	#else
903	if ((de = (int)(d0 >> Exp_shift)))
904	z |= Exp_msk1;
905	#endif
906	#ifdef Pack_32
907	if ((y = d1)) {
908	if ((k = lo0bits(y: &y))) {
909	x[0] = y | (z << (32 - k));
910	z >>= k;
911	} else
912	x[0] = y;
913	if (z) {
914	b.resize(s: 2);
915	x[1] = z;
916	}
917
918	#ifndef Sudden_Underflow
919	i = b.size();
920	#endif
921	} else {
922	k = lo0bits(y: &z);
923	x[0] = z;
924	#ifndef Sudden_Underflow
925	i = 1;
926	#endif
927	b.resize(s: 1);
928	k += 32;
929	}
930	#else
931	if ((y = d1)) {
932	if ((k = lo0bits(&y))) {
933	if (k >= 16) {
934	x[0] = y | z << 32 - k & 0xffff;
935	x[1] = z >> k - 16 & 0xffff;
936	x[2] = z >> k;
937	i = 2;
938	} else {
939	x[0] = y & 0xffff;
940	x[1] = y >> 16 | z << 16 - k & 0xffff;
941	x[2] = z >> k & 0xffff;
942	x[3] = z >> k + 16;
943	i = 3;
944	}
945	} else {
946	x[0] = y & 0xffff;
947	x[1] = y >> 16;
948	x[2] = z & 0xffff;
949	x[3] = z >> 16;
950	i = 3;
951	}
952	} else {
953	k = lo0bits(&z);
954	if (k >= 16) {
955	x[0] = z;
956	i = 0;
957	} else {
958	x[0] = z & 0xffff;
959	x[1] = z >> 16;
960	i = 1;
961	}
962	k += 32;
963	} while (!x[i])
964	--i;
965	b->resize(i + 1);
966	#endif
967	#ifndef Sudden_Underflow
968	if (de) {
969	#endif
970	*e = de - Bias - (P - 1) + k;
971	*bits = P - k;
972	#ifndef Sudden_Underflow
973	} else {
974	*e = de - Bias - (P - 1) + 1 + k;
975	#ifdef Pack_32
976	*bits = (32 * i) - hi0bits(x: x[i - 1]);
977	#else
978	*bits = (i + 2) * 16 - hi0bits(x[i]);
979	#endif
980	}
981	#endif
982	}
983	#undef d0
984	#undef d1
985
986	static double ratio(const BigInt& a, const BigInt& b)
987	{
988	U da, db;
989	int k, ka, kb;
990
991	dval(&da) = b2d(a, e: &ka);
992	dval(&db) = b2d(a: b, e: &kb);
993	#ifdef Pack_32
994	k = ka - kb + 32 * (a.size() - b.size());
995	#else
996	k = ka - kb + 16 * (a.size() - b.size());
997	#endif
998	if (k > 0)
999	word0(&da) += k * Exp_msk1;
1000	else {
1001	k = -k;
1002	word0(&db) += k * Exp_msk1;
1003	}
1004	return dval(&da) / dval(&db);
1005	}
1006
1007	static const double tens[] = {
1008	1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1009	1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1010	1e20, 1e21, 1e22
1011	};
1012
1013	static const double bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1014	static const double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
1015	#ifdef Avoid_Underflow
1016	9007199254740992. * 9007199254740992.e-256
1017	/* = 2^106 * 1e-53 */
1018	#else
1019	1e-256
1020	#endif
1021	};
1022
1023	/* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
1024	/* flag unnecessarily. It leads to a song and dance at the end of strtod. */
1025	#define Scale_Bit 0x10
1026	#define n_bigtens 5
1027
1028	#if defined(INFNAN_CHECK)
1029
1030	#ifndef NAN_WORD0
1031	#define NAN_WORD0 0x7ff80000
1032	#endif
1033
1034	#ifndef NAN_WORD1
1035	#define NAN_WORD1 0
1036	#endif
1037
1038	static int match(const char** sp, const char* t)
1039	{
1040	int c, d;
1041	const char* s = *sp;
1042
1043	while ((d = *t++)) {
1044	if ((c = *++s) >= 'A' && c <= 'Z')
1045	c += 'a' - 'A';
1046	if (c != d)
1047	return 0;
1048	}
1049	*sp = s + 1;
1050	return 1;
1051	}
1052
1053	#ifndef No_Hex_NaN
1054	static void hexnan(U* rvp, const char** sp)
1055	{
1056	uint32_t c, x[2];
1057	const char* s;
1058	int havedig, udx0, xshift;
1059
1060	x[0] = x[1] = 0;
1061	havedig = xshift = 0;
1062	udx0 = 1;
1063	s = *sp;
1064	while ((c = *(const unsigned char*)++s)) {
1065	if (c >= '0' && c <= '9')
1066	c -= '0';
1067	else if (c >= 'a' && c <= 'f')
1068	c += 10 - 'a';
1069	else if (c >= 'A' && c <= 'F')
1070	c += 10 - 'A';
1071	else if (c <= ' ') {
1072	if (udx0 && havedig) {
1073	udx0 = 0;
1074	xshift = 1;
1075	}
1076	continue;
1077	} else if (/*(*/ c == ')' && havedig) {
1078	*sp = s + 1;
1079	break;
1080	} else
1081	return; /* invalid form: don't change *sp */
1082	havedig = 1;
1083	if (xshift) {
1084	xshift = 0;
1085	x[0] = x[1];
1086	x[1] = 0;
1087	}
1088	if (udx0)
1089	x[0] = (x[0] << 4) | (x[1] >> 28);
1090	x[1] = (x[1] << 4) | c;
1091	}
1092	if ((x[0] &= 0xfffff) || x[1]) {
1093	word0(rvp) = Exp_mask | x[0];
1094	word1(rvp) = x[1];
1095	}
1096	}
1097	#endif /*No_Hex_NaN*/
1098	#endif /* INFNAN_CHECK */
1099
1100	double strtod(const char* s00, char** se)
1101	{
1102	#ifdef Avoid_Underflow
1103	int scale;
1104	#endif
1105	int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
1106	e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
1107	const char *s, *s0, *s1;
1108	double aadj, aadj1;
1109	U aadj2, adj, rv, rv0;
1110	int32_t L;
1111	uint32_t y, z;
1112	BigInt bb, bb1, bd, bd0, bs, delta;
1113	#ifdef SET_INEXACT
1114	int inexact, oldinexact;
1115	#endif
1116
1117	sign = nz0 = nz = 0;
1118	dval(&rv) = 0;
1119	for (s = s00; ; s++)
1120	switch (*s) {
1121	case '-':
1122	sign = 1;
1123	/* no break */
1124	case '+':
1125	if (*++s)
1126	goto break2;
1127	/* no break */
1128	case 0:
1129	goto ret0;
1130	case '\t':
1131	case '\n':
1132	case '\v':
1133	case '\f':
1134	case '\r':
1135	case ' ':
1136	continue;
1137	default:
1138	goto break2;
1139	}
1140	break2:
1141	if (*s == '0') {
1142	nz0 = 1;
1143	while (*++s == '0') { }
1144	if (!*s)
1145	goto ret;
1146	}
1147	s0 = s;
1148	y = z = 0;
1149	for (nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
1150	if (nd < 9)
1151	y = (10 * y) + c - '0';
1152	else if (nd < 16)
1153	z = (10 * z) + c - '0';
1154	nd0 = nd;
1155	if (c == '.') {
1156	c = *++s;
1157	if (!nd) {
1158	for (; c == '0'; c = *++s)
1159	nz++;
1160	if (c > '0' && c <= '9') {
1161	s0 = s;
1162	nf += nz;
1163	nz = 0;
1164	goto have_dig;
1165	}
1166	goto dig_done;
1167	}
1168	for (; c >= '0' && c <= '9'; c = *++s) {
1169	have_dig:
1170	nz++;
1171	if (c -= '0') {
1172	nf += nz;
1173	for (i = 1; i < nz; i++)
1174	if (nd++ < 9)
1175	y *= 10;
1176	else if (nd <= DBL_DIG + 1)
1177	z *= 10;
1178	if (nd++ < 9)
1179	y = (10 * y) + c;
1180	else if (nd <= DBL_DIG + 1)
1181	z = (10 * z) + c;
1182	nz = 0;
1183	}
1184	}
1185	}
1186	dig_done:
1187	e = 0;
1188	if (c == 'e' || c == 'E') {
1189	if (!nd && !nz && !nz0) {
1190	goto ret0;
1191	}
1192	s00 = s;
1193	esign = 0;
1194	switch (c = *++s) {
1195	case '-':
1196	esign = 1;
1197	case '+':
1198	c = *++s;
1199	}
1200	if (c >= '0' && c <= '9') {
1201	while (c == '0')
1202	c = *++s;
1203	if (c > '0' && c <= '9') {
1204	L = c - '0';
1205	s1 = s;
1206	while ((c = *++s) >= '0' && c <= '9')
1207	L = (10 * L) + c - '0';
1208	if (s - s1 > 8 || L > 19999)
1209	/* Avoid confusion from exponents
1210	* so large that e might overflow.
1211	*/
1212	e = 19999; /* safe for 16 bit ints */
1213	else
1214	e = (int)L;
1215	if (esign)
1216	e = -e;
1217	} else
1218	e = 0;
1219	} else
1220	s = s00;
1221	}
1222	if (!nd) {
1223	if (!nz && !nz0) {
1224	#ifdef INFNAN_CHECK
1225	/* Check for Nan and Infinity */
1226	switch(c) {
1227	case 'i':
1228	case 'I':
1229	if (match(sp: &s,t: "nf")) {
1230	--s;
1231	if (!match(sp: &s,t: "inity"))
1232	++s;
1233	word0(&rv) = 0x7ff00000;
1234	word1(&rv) = 0;
1235	goto ret;
1236	}
1237	break;
1238	case 'n':
1239	case 'N':
1240	if (match(sp: &s, t: "an")) {
1241	word0(&rv) = NAN_WORD0;
1242	word1(&rv) = NAN_WORD1;
1243	#ifndef No_Hex_NaN
1244	if (*s == '(') /*)*/
1245	hexnan(rvp: &rv, sp: &s);
1246	#endif
1247	goto ret;
1248	}
1249	}
1250	#endif /* INFNAN_CHECK */
1251	ret0:
1252	s = s00;
1253	sign = 0;
1254	}
1255	goto ret;
1256	}
1257	e1 = e -= nf;
1258
1259	/* Now we have nd0 digits, starting at s0, followed by a
1260	* decimal point, followed by nd-nd0 digits. The number we're
1261	* after is the integer represented by those digits times
1262	* 10**e */
1263
1264	if (!nd0)
1265	nd0 = nd;
1266	k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
1267	dval(&rv) = y;
1268	if (k > 9) {
1269	#ifdef SET_INEXACT
1270	if (k > DBL_DIG)
1271	oldinexact = get_inexact();
1272	#endif
1273	dval(&rv) = tens[k - 9] * dval(&rv) + z;
1274	}
1275	if (nd <= DBL_DIG && Flt_Rounds == 1) {
1276	if (!e)
1277	goto ret;
1278	if (e > 0) {
1279	if (e <= Ten_pmax) {
1280	/* rv = */ rounded_product(dval(&rv), tens[e]);
1281	goto ret;
1282	}
1283	i = DBL_DIG - nd;
1284	if (e <= Ten_pmax + i) {
1285	/* A fancier test would sometimes let us do
1286	* this for larger i values.
1287	*/
1288	e -= i;
1289	dval(&rv) *= tens[i];
1290	/* rv = */ rounded_product(dval(&rv), tens[e]);
1291	goto ret;
1292	}
1293	}
1294	#ifndef Inaccurate_Divide
1295	else if (e >= -Ten_pmax) {
1296	/* rv = */ rounded_quotient(dval(&rv), tens[-e]);
1297	goto ret;
1298	}
1299	#endif
1300	}
1301	e1 += nd - k;
1302
1303	#ifdef SET_INEXACT
1304	inexact = 1;
1305	if (k <= DBL_DIG)
1306	oldinexact = get_inexact();
1307	#endif
1308	#ifdef Avoid_Underflow
1309	scale = 0;
1310	#endif
1311
1312	/* Get starting approximation = rv * 10**e1 */
1313
1314	if (e1 > 0) {
1315	if ((i = e1 & 15))
1316	dval(&rv) *= tens[i];
1317	if (e1 &= ~15) {
1318	if (e1 > DBL_MAX_10_EXP) {
1319	ovfl:
1320	#ifndef NO_ERRNO
1321	errno = ERANGE;
1322	#endif
1323	/* Can't trust HUGE_VAL */
1324	word0(&rv) = Exp_mask;
1325	word1(&rv) = 0;
1326	#ifdef SET_INEXACT
1327	/* set overflow bit */
1328	dval(&rv0) = 1e300;
1329	dval(&rv0) *= dval(&rv0);
1330	#endif
1331	goto ret;
1332	}
1333	e1 >>= 4;
1334	for (j = 0; e1 > 1; j++, e1 >>= 1)
1335	if (e1 & 1)
1336	dval(&rv) *= bigtens[j];
1337	/* The last multiplication could overflow. */
1338	word0(&rv) -= P * Exp_msk1;
1339	dval(&rv) *= bigtens[j];
1340	if ((z = word0(&rv) & Exp_mask) > Exp_msk1 * (DBL_MAX_EXP + Bias - P))
1341	goto ovfl;
1342	if (z > Exp_msk1 * (DBL_MAX_EXP + Bias - 1 - P)) {
1343	/* set to largest number */
1344	/* (Can't trust DBL_MAX) */
1345	word0(&rv) = Big0;
1346	word1(&rv) = Big1;
1347	} else
1348	word0(&rv) += P * Exp_msk1;
1349	}
1350	} else if (e1 < 0) {
1351	e1 = -e1;
1352	if ((i = e1 & 15))
1353	dval(&rv) /= tens[i];
1354	if (e1 >>= 4) {
1355	if (e1 >= 1 << n_bigtens)
1356	goto undfl;
1357	#ifdef Avoid_Underflow
1358	if (e1 & Scale_Bit)
1359	scale = 2 * P;
1360	for (j = 0; e1 > 0; j++, e1 >>= 1)
1361	if (e1 & 1)
1362	dval(&rv) *= tinytens[j];
1363	if (scale && (j = (2 * P) + 1 - ((word0(&rv) & Exp_mask) >> Exp_shift)) > 0) {
1364	/* scaled rv is denormal; zap j low bits */
1365	if (j >= 32) {
1366	word1(&rv) = 0;
1367	if (j >= 53)
1368	word0(&rv) = (P + 2) * Exp_msk1;
1369	else
1370	word0(&rv) &= 0xffffffff << (j - 32);
1371	} else
1372	word1(&rv) &= 0xffffffff << j;
1373	}
1374	#else
1375	for (j = 0; e1 > 1; j++, e1 >>= 1)
1376	if (e1 & 1)
1377	dval(&rv) *= tinytens[j];
1378	/* The last multiplication could underflow. */
1379	dval(&rv0) = dval(&rv);
1380	dval(&rv) *= tinytens[j];
1381	if (!dval(&rv)) {
1382	dval(&rv) = 2. * dval(&rv0);
1383	dval(&rv) *= tinytens[j];
1384	#endif
1385	if (!dval(&rv)) {
1386	undfl:
1387	dval(&rv) = 0.;
1388	#ifndef NO_ERRNO
1389	errno = ERANGE;
1390	#endif
1391	goto ret;
1392	}
1393	#ifndef Avoid_Underflow
1394	word0(&rv) = Tiny0;
1395	word1(&rv) = Tiny1;
1396	/* The refinement below will clean
1397	* this approximation up.
1398	*/
1399	}
1400	#endif
1401	}
1402	}
1403
1404	/* Now the hard part -- adjusting rv to the correct value.*/
1405
1406	/* Put digits into bd: true value = bd * 10^e */
1407
1408	s2b(b&: bd0, s: s0, nd0, nd, y9: y);
1409
1410	for (;;) {
1411	bd = bd0;
1412	d2b(b&: bb, d: &rv, e: &bbe, bits: &bbbits); /* rv = bb * 2^bbe */
1413	i2b(b&: bs, i: 1);
1414
1415	if (e >= 0) {
1416	bb2 = bb5 = 0;
1417	bd2 = bd5 = e;
1418	} else {
1419	bb2 = bb5 = -e;
1420	bd2 = bd5 = 0;
1421	}
1422	if (bbe >= 0)
1423	bb2 += bbe;
1424	else
1425	bd2 -= bbe;
1426	bs2 = bb2;
1427	#ifdef Avoid_Underflow
1428	j = bbe - scale;
1429	i = j + bbbits - 1; /* logb(rv) */
1430	if (i < Emin) /* denormal */
1431	j += P - Emin;
1432	else
1433	j = P + 1 - bbbits;
1434	#else /*Avoid_Underflow*/
1435	#ifdef Sudden_Underflow
1436	j = P + 1 - bbbits;
1437	#else /*Sudden_Underflow*/
1438	j = bbe;
1439	i = j + bbbits - 1; /* logb(rv) */
1440	if (i < Emin) /* denormal */
1441	j += P - Emin;
1442	else
1443	j = P + 1 - bbbits;
1444	#endif /*Sudden_Underflow*/
1445	#endif /*Avoid_Underflow*/
1446	bb2 += j;
1447	bd2 += j;
1448	#ifdef Avoid_Underflow
1449	bd2 += scale;
1450	#endif
1451	i = bb2 < bd2 ? bb2 : bd2;
1452	if (i > bs2)
1453	i = bs2;
1454	if (i > 0) {
1455	bb2 -= i;
1456	bd2 -= i;
1457	bs2 -= i;
1458	}
1459	if (bb5 > 0) {
1460	pow5mult(b&: bs, k: bb5);
1461	mult(aRef&: bb, bRef: bs);
1462	}
1463	if (bb2 > 0)
1464	lshift(b&: bb, k: bb2);
1465	if (bd5 > 0)
1466	pow5mult(b&: bd, k: bd5);
1467	if (bd2 > 0)
1468	lshift(b&: bd, k: bd2);
1469	if (bs2 > 0)
1470	lshift(b&: bs, k: bs2);
1471	diff(c&: delta, aRef: bb, bRef: bd);
1472	dsign = delta.sign;
1473	delta.sign = 0;
1474	i = cmp(a: delta, b: bs);
1475
1476	if (i < 0) {
1477	/* Error is less than half an ulp -- check for
1478	* special case of mantissa a power of two.
1479	*/
1480	if (dsign || word1(&rv) || word0(&rv) & Bndry_mask
1481	#ifdef Avoid_Underflow
1482	|| (word0(&rv) & Exp_mask) <= (2 * P + 1) * Exp_msk1
1483	#else
1484	|| (word0(&rv) & Exp_mask) <= Exp_msk1
1485	#endif
1486	) {
1487	#ifdef SET_INEXACT
1488	if (!delta->words()[0] && delta->size() <= 1)
1489	inexact = 0;
1490	#endif
1491	break;
1492	}
1493	if (!delta.words()[0] && delta.size() <= 1) {
1494	/* exact result */
1495	#ifdef SET_INEXACT
1496	inexact = 0;
1497	#endif
1498	break;
1499	}
1500	lshift(b&: delta, Log2P);
1501	if (cmp(a: delta, b: bs) > 0)
1502	goto drop_down;
1503	break;
1504	}
1505	if (i == 0) {
1506	/* exactly half-way between */
1507	if (dsign) {
1508	if ((word0(&rv) & Bndry_mask1) == Bndry_mask1
1509	&& word1(&rv) == (
1510	#ifdef Avoid_Underflow
1511	(scale && (y = word0(&rv) & Exp_mask) <= 2 * P * Exp_msk1)
1512	? (0xffffffff & (0xffffffff << (2 * P + 1 - (y >> Exp_shift)))) :
1513	#endif
1514	0xffffffff)) {
1515	/*boundary case -- increment exponent*/
1516	word0(&rv) = (word0(&rv) & Exp_mask) + Exp_msk1;
1517	word1(&rv) = 0;
1518	#ifdef Avoid_Underflow
1519	dsign = 0;
1520	#endif
1521	break;
1522	}
1523	} else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) {
1524	drop_down:
1525	/* boundary case -- decrement exponent */
1526	#ifdef Sudden_Underflow /*{{*/
1527	L = word0(&rv) & Exp_mask;
1528	#ifdef Avoid_Underflow
1529	if (L <= (scale ? (2 * P + 1) * Exp_msk1 : Exp_msk1))
1530	#else
1531	if (L <= Exp_msk1)
1532	#endif /*Avoid_Underflow*/
1533	goto undfl;
1534	L -= Exp_msk1;
1535	#else /*Sudden_Underflow}{*/
1536	#ifdef Avoid_Underflow
1537	if (scale) {
1538	L = word0(&rv) & Exp_mask;
1539	if (L <= (2 * P + 1) * Exp_msk1) {
1540	if (L > (P + 2) * Exp_msk1)
1541	/* round even ==> */
1542	/* accept rv */
1543	break;
1544	/* rv = smallest denormal */
1545	goto undfl;
1546	}
1547	}
1548	#endif /*Avoid_Underflow*/
1549	L = (word0(&rv) & Exp_mask) - Exp_msk1;
1550	#endif /*Sudden_Underflow}}*/
1551	word0(&rv) = L | Bndry_mask1;
1552	word1(&rv) = 0xffffffff;
1553	break;
1554	}
1555	if (!(word1(&rv) & LSB))
1556	break;
1557	if (dsign)
1558	dval(&rv) += ulp(x: &rv);
1559	else {
1560	dval(&rv) -= ulp(x: &rv);
1561	#ifndef Sudden_Underflow
1562	if (!dval(&rv))
1563	goto undfl;
1564	#endif
1565	}
1566	#ifdef Avoid_Underflow
1567	dsign = 1 - dsign;
1568	#endif
1569	break;
1570	}
1571	if ((aadj = ratio(a: delta, b: bs)) <= 2.) {
1572	if (dsign)
1573	aadj = aadj1 = 1.;
1574	else if (word1(&rv) || word0(&rv) & Bndry_mask) {
1575	#ifndef Sudden_Underflow
1576	if (word1(&rv) == Tiny1 && !word0(&rv))
1577	goto undfl;
1578	#endif
1579	aadj = 1.;
1580	aadj1 = -1.;
1581	} else {
1582	/* special case -- power of FLT_RADIX to be */
1583	/* rounded down... */
1584
1585	if (aadj < 2. / FLT_RADIX)
1586	aadj = 1. / FLT_RADIX;
1587	else
1588	aadj *= 0.5;
1589	aadj1 = -aadj;
1590	}
1591	} else {
1592	aadj *= 0.5;
1593	aadj1 = dsign ? aadj : -aadj;
1594	#ifdef Check_FLT_ROUNDS
1595	switch (Rounding) {
1596	case 2: /* towards +infinity */
1597	aadj1 -= 0.5;
1598	break;
1599	case 0: /* towards 0 */
1600	case 3: /* towards -infinity */
1601	aadj1 += 0.5;
1602	}
1603	#else
1604	if (Flt_Rounds == 0)
1605	aadj1 += 0.5;
1606	#endif /*Check_FLT_ROUNDS*/
1607	}
1608	y = word0(&rv) & Exp_mask;
1609
1610	/* Check for overflow */
1611
1612	if (y == Exp_msk1 * (DBL_MAX_EXP + Bias - 1)) {
1613	dval(&rv0) = dval(&rv);
1614	word0(&rv) -= P * Exp_msk1;
1615	adj.d = aadj1 * ulp(x: &rv);
1616	dval(&rv) += adj.d;
1617	if ((word0(&rv) & Exp_mask) >= Exp_msk1 * (DBL_MAX_EXP + Bias - P)) {
1618	if (word0(&rv0) == Big0 && word1(&rv0) == Big1)
1619	goto ovfl;
1620	word0(&rv) = Big0;
1621	word1(&rv) = Big1;
1622	goto cont;
1623	} else
1624	word0(&rv) += P * Exp_msk1;
1625	} else {
1626	#ifdef Avoid_Underflow
1627	if (scale && y <= 2 * P * Exp_msk1) {
1628	if (aadj <= 0x7fffffff) {
1629	if ((z = (uint32_t)aadj) <= 0)
1630	z = 1;
1631	aadj = z;
1632	aadj1 = dsign ? aadj : -aadj;
1633	}
1634	dval(&aadj2) = aadj1;
1635	word0(&aadj2) += (2 * P + 1) * Exp_msk1 - y;
1636	aadj1 = dval(&aadj2);
1637	}
1638	adj.d = aadj1 * ulp(x: &rv);
1639	dval(&rv) += adj.d;
1640	#else
1641	#ifdef Sudden_Underflow
1642	if ((word0(&rv) & Exp_mask) <= P * Exp_msk1) {
1643	dval(&rv0) = dval(&rv);
1644	word0(&rv) += P * Exp_msk1;
1645	adj.d = aadj1 * ulp(&rv);
1646	dval(&rv) += adj.d;
1647	if ((word0(&rv) & Exp_mask) <= P * Exp_msk1)
1648	{
1649	if (word0(&rv0) == Tiny0 && word1(&rv0) == Tiny1)
1650	goto undfl;
1651	word0(&rv) = Tiny0;
1652	word1(&rv) = Tiny1;
1653	goto cont;
1654	}
1655	else
1656	word0(&rv) -= P * Exp_msk1;
1657	} else {
1658	adj.d = aadj1 * ulp(&rv);
1659	dval(&rv) += adj.d;
1660	}
1661	#else /*Sudden_Underflow*/
1662	/* Compute adj so that the IEEE rounding rules will
1663	* correctly round rv + adj in some half-way cases.
1664	* If rv * ulp(rv) is denormalized (i.e.,
1665	* y <= (P - 1) * Exp_msk1), we must adjust aadj to avoid
1666	* trouble from bits lost to denormalization;
1667	* example: 1.2e-307 .
1668	*/
1669	if (y <= (P - 1) * Exp_msk1 && aadj > 1.) {
1670	aadj1 = (double)(int)(aadj + 0.5);
1671	if (!dsign)
1672	aadj1 = -aadj1;
1673	}
1674	adj.d = aadj1 * ulp(&rv);
1675	dval(&rv) += adj.d;
1676	#endif /*Sudden_Underflow*/
1677	#endif /*Avoid_Underflow*/
1678	}
1679	z = word0(&rv) & Exp_mask;
1680	#ifndef SET_INEXACT
1681	#ifdef Avoid_Underflow
1682	if (!scale)
1683	#endif
1684	if (y == z) {
1685	/* Can we stop now? */
1686	L = (int32_t)aadj;
1687	aadj -= L;
1688	/* The tolerances below are conservative. */
1689	if (dsign || word1(&rv) || word0(&rv) & Bndry_mask) {
1690	if (aadj < .4999999 || aadj > .5000001)
1691	break;
1692	} else if (aadj < .4999999 / FLT_RADIX)
1693	break;
1694	}
1695	#endif
1696	cont:
1697	;
1698	}
1699	#ifdef SET_INEXACT
1700	if (inexact) {
1701	if (!oldinexact) {
1702	word0(&rv0) = Exp_1 + (70 << Exp_shift);
1703	word1(&rv0) = 0;
1704	dval(&rv0) += 1.;
1705	}
1706	} else if (!oldinexact)
1707	clear_inexact();
1708	#endif
1709	#ifdef Avoid_Underflow
1710	if (scale) {
1711	word0(&rv0) = Exp_1 - 2 * P * Exp_msk1;
1712	word1(&rv0) = 0;
1713	dval(&rv) *= dval(&rv0);
1714	#ifndef NO_ERRNO
1715	/* try to avoid the bug of testing an 8087 register value */
1716	if (word0(&rv) == 0 && word1(&rv) == 0)
1717	errno = ERANGE;
1718	#endif
1719	}
1720	#endif /* Avoid_Underflow */
1721	#ifdef SET_INEXACT
1722	if (inexact && !(word0(&rv) & Exp_mask)) {
1723	/* set underflow bit */
1724	dval(&rv0) = 1e-300;
1725	dval(&rv0) *= dval(&rv0);
1726	}
1727	#endif
1728	ret:
1729	if (se)
1730	*se = const_cast<char*>(s);
1731	return sign ? -dval(&rv) : dval(&rv);
1732	}
1733
1734	static ALWAYS_INLINE int quorem(BigInt& b, BigInt& S)
1735	{
1736	size_t n;
1737	uint32_t *bx, *bxe, q, *sx, *sxe;
1738	#ifdef USE_LONG_LONG
1739	unsigned long long borrow, carry, y, ys;
1740	#else
1741	uint32_t borrow, carry, y, ys;
1742	#ifdef Pack_32
1743	uint32_t si, z, zs;
1744	#endif
1745	#endif
1746	ASSERT(b.size() <= 1 || b.words()[b.size() - 1]);
1747	ASSERT(S.size() <= 1 || S.words()[S.size() - 1]);
1748
1749	n = S.size();
1750	ASSERT_WITH_MESSAGE(b.size() <= n, "oversize b in quorem");
1751	if (b.size() < n)
1752	return 0;
1753	sx = S.words();
1754	sxe = sx + --n;
1755	bx = b.words();
1756	bxe = bx + n;
1757	q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
1758	ASSERT_WITH_MESSAGE(q <= 9, "oversized quotient in quorem");
1759	if (q) {
1760	borrow = 0;
1761	carry = 0;
1762	do {
1763	#ifdef USE_LONG_LONG
1764	ys = *sx++ * (unsigned long long)q + carry;
1765	carry = ys >> 32;
1766	y = *bx - (ys & 0xffffffffUL) - borrow;
1767	borrow = y >> 32 & (uint32_t)1;
1768	*bx++ = (uint32_t)y & 0xffffffffUL;
1769	#else
1770	#ifdef Pack_32
1771	si = *sx++;
1772	ys = (si & 0xffff) * q + carry;
1773	zs = (si >> 16) * q + (ys >> 16);
1774	carry = zs >> 16;
1775	y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
1776	borrow = (y & 0x10000) >> 16;
1777	z = (*bx >> 16) - (zs & 0xffff) - borrow;
1778	borrow = (z & 0x10000) >> 16;
1779	Storeinc(bx, z, y);
1780	#else
1781	ys = *sx++ * q + carry;
1782	carry = ys >> 16;
1783	y = *bx - (ys & 0xffff) - borrow;
1784	borrow = (y & 0x10000) >> 16;
1785	*bx++ = y & 0xffff;
1786	#endif
1787	#endif
1788	} while (sx <= sxe);
1789	if (!*bxe) {
1790	bx = b.words();
1791	while (--bxe > bx && !*bxe)
1792	--n;
1793	b.resize(s: n);
1794	}
1795	}
1796	if (cmp(a: b, b: S) >= 0) {
1797	q++;
1798	borrow = 0;
1799	carry = 0;
1800	bx = b.words();
1801	sx = S.words();
1802	do {
1803	#ifdef USE_LONG_LONG
1804	ys = *sx++ + carry;
1805	carry = ys >> 32;
1806	y = *bx - (ys & 0xffffffffUL) - borrow;
1807	borrow = y >> 32 & (uint32_t)1;
1808	*bx++ = (uint32_t)y & 0xffffffffUL;
1809	#else
1810	#ifdef Pack_32
1811	si = *sx++;
1812	ys = (si & 0xffff) + carry;
1813	zs = (si >> 16) + (ys >> 16);
1814	carry = zs >> 16;
1815	y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
1816	borrow = (y & 0x10000) >> 16;
1817	z = (*bx >> 16) - (zs & 0xffff) - borrow;
1818	borrow = (z & 0x10000) >> 16;
1819	Storeinc(bx, z, y);
1820	#else
1821	ys = *sx++ + carry;
1822	carry = ys >> 16;
1823	y = *bx - (ys & 0xffff) - borrow;
1824	borrow = (y & 0x10000) >> 16;
1825	*bx++ = y & 0xffff;
1826	#endif
1827	#endif
1828	} while (sx <= sxe);
1829	bx = b.words();
1830	bxe = bx + n;
1831	if (!*bxe) {
1832	while (--bxe > bx && !*bxe)
1833	--n;
1834	b.resize(s: n);
1835	}
1836	}
1837	return q;
1838	}
1839
1840	/* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
1841	*
1842	* Inspired by "How to Print Floating-Point Numbers Accurately" by
1843	* Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
1844	*
1845	* Modifications:
1846	* 1. Rather than iterating, we use a simple numeric overestimate
1847	* to determine k = floor(log10(d)). We scale relevant
1848	* quantities using O(log2(k)) rather than O(k) multiplications.
1849	* 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
1850	* try to generate digits strictly left to right. Instead, we
1851	* compute with fewer bits and propagate the carry if necessary
1852	* when rounding the final digit up. This is often faster.
1853	* 3. Under the assumption that input will be rounded nearest,
1854	* mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
1855	* That is, we allow equality in stopping tests when the
1856	* round-nearest rule will give the same floating-point value
1857	* as would satisfaction of the stopping test with strict
1858	* inequality.
1859	* 4. We remove common factors of powers of 2 from relevant
1860	* quantities.
1861	* 5. When converting floating-point integers less than 1e16,
1862	* we use floating-point arithmetic rather than resorting
1863	* to multiple-precision integers.
1864	* 6. When asked to produce fewer than 15 digits, we first try
1865	* to get by with floating-point arithmetic; we resort to
1866	* multiple-precision integer arithmetic only if we cannot
1867	* guarantee that the floating-point calculation has given
1868	* the correctly rounded result. For k requested digits and
1869	* "uniformly" distributed input, the probability is
1870	* something like 10^(k-15) that we must resort to the int32_t
1871	* calculation.
1872	*/
1873
1874	void dtoa(DtoaBuffer result, double dd, int ndigits, int* decpt, int* sign, char** rve)
1875	{
1876	/*
1877	Arguments ndigits, decpt, sign are similar to those
1878	of ecvt and fcvt; trailing zeros are suppressed from
1879	the returned string. If not null, *rve is set to point
1880	to the end of the return value. If d is +-Infinity or NaN,
1881	then *decpt is set to 9999.
1882
1883	*/
1884
1885	int bbits, b2, b5, be, dig, i, ieps, ilim = 0, ilim0, ilim1 = 0,
1886	j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
1887	spec_case, try_quick;
1888	int32_t L;
1889	#ifndef Sudden_Underflow
1890	int denorm;
1891	uint32_t x;
1892	#endif
1893	BigInt b, b1, delta, mlo, mhi, S;
1894	U d2, eps, u;
1895	double ds;
1896	char *s, *s0;
1897	#ifdef SET_INEXACT
1898	int inexact, oldinexact;
1899	#endif
1900
1901	u.d = dd;
1902	if (word0(&u) & Sign_bit) {
1903	/* set sign for everything, including 0's and NaNs */
1904	*sign = 1;
1905	word0(&u) &= ~Sign_bit; /* clear sign bit */
1906	} else
1907	*sign = 0;
1908
1909	if ((word0(&u) & Exp_mask) == Exp_mask)
1910	{
1911	/* Infinity or NaN */
1912	*decpt = 9999;
1913	if (!word1(&u) && !(word0(&u) & 0xfffff)) {
1914	strcpy(dest: result, src: "Infinity");
1915	if (rve)
1916	*rve = result + 8;
1917	} else {
1918	strcpy(dest: result, src: "NaN");
1919	if (rve)
1920	*rve = result + 3;
1921	}
1922	return;
1923	}
1924	if (!dval(&u)) {
1925	*decpt = 1;
1926	result[0] = '0';
1927	result[1] = '\0';
1928	if (rve)
1929	*rve = result + 1;
1930	return;
1931	}
1932
1933	#ifdef SET_INEXACT
1934	try_quick = oldinexact = get_inexact();
1935	inexact = 1;
1936	#endif
1937
1938	d2b(b, d: &u, e: &be, bits: &bbits);
1939	#ifdef Sudden_Underflow
1940	i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask >> Exp_shift1));
1941	#else
1942	if ((i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask >> Exp_shift1)))) {
1943	#endif
1944	dval(&d2) = dval(&u);
1945	word0(&d2) &= Frac_mask1;
1946	word0(&d2) |= Exp_11;
1947
1948	/* log(x) ~=~ log(1.5) + (x-1.5)/1.5
1949	* log10(x) = log(x) / log(10)
1950	* ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
1951	* log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
1952	*
1953	* This suggests computing an approximation k to log10(d) by
1954	*
1955	* k = (i - Bias)*0.301029995663981
1956	* + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
1957	*
1958	* We want k to be too large rather than too small.
1959	* The error in the first-order Taylor series approximation
1960	* is in our favor, so we just round up the constant enough
1961	* to compensate for any error in the multiplication of
1962	* (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
1963	* and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
1964	* adding 1e-13 to the constant term more than suffices.
1965	* Hence we adjust the constant term to 0.1760912590558.
1966	* (We could get a more accurate k by invoking log10,
1967	* but this is probably not worthwhile.)
1968	*/
1969
1970	i -= Bias;
1971	#ifndef Sudden_Underflow
1972	denorm = 0;
1973	} else {
1974	/* d is denormalized */
1975
1976	i = bbits + be + (Bias + (P - 1) - 1);
1977	x = (i > 32) ? (word0(&u) << (64 - i)) | (word1(&u) >> (i - 32))
1978	: word1(&u) << (32 - i);
1979	dval(&d2) = x;
1980	word0(&d2) -= 31 * Exp_msk1; /* adjust exponent */
1981	i -= (Bias + (P - 1) - 1) + 1;
1982	denorm = 1;
1983	}
1984	#endif
1985	ds = (dval(&d2) - 1.5) * 0.289529654602168 + 0.1760912590558 + (i * 0.301029995663981);
1986	k = (int)ds;
1987	if (ds < 0. && ds != k)
1988	k--; /* want k = floor(ds) */
1989	k_check = 1;
1990	if (k >= 0 && k <= Ten_pmax) {
1991	if (dval(&u) < tens[k])
1992	k--;
1993	k_check = 0;
1994	}
1995	j = bbits - i - 1;
1996	if (j >= 0) {
1997	b2 = 0;
1998	s2 = j;
1999	} else {
2000	b2 = -j;
2001	s2 = 0;
2002	}
2003	if (k >= 0) {
2004	b5 = 0;
2005	s5 = k;
2006	s2 += k;
2007	} else {
2008	b2 -= k;
2009	b5 = -k;
2010	s5 = 0;
2011	}
2012
2013	#ifndef SET_INEXACT
2014	#ifdef Check_FLT_ROUNDS
2015	try_quick = Rounding == 1;
2016	#else
2017	try_quick = 1;
2018	#endif
2019	#endif /*SET_INEXACT*/
2020
2021	leftright = 1;
2022	ilim = ilim1 = -1;
2023	i = 18;
2024	ndigits = 0;
2025	s = s0 = result;
2026
2027	if (ilim >= 0 && ilim <= Quick_max && try_quick) {
2028
2029	/* Try to get by with floating-point arithmetic. */
2030
2031	i = 0;
2032	dval(&d2) = dval(&u);
2033	k0 = k;
2034	ilim0 = ilim;
2035	ieps = 2; /* conservative */
2036	if (k > 0) {
2037	ds = tens[k & 0xf];
2038	j = k >> 4;
2039	if (j & Bletch) {
2040	/* prevent overflows */
2041	j &= Bletch - 1;
2042	dval(&u) /= bigtens[n_bigtens - 1];
2043	ieps++;
2044	}
2045	for (; j; j >>= 1, i++) {
2046	if (j & 1) {
2047	ieps++;
2048	ds *= bigtens[i];
2049	}
2050	}
2051	dval(&u) /= ds;
2052	} else if ((j1 = -k)) {
2053	dval(&u) *= tens[j1 & 0xf];
2054	for (j = j1 >> 4; j; j >>= 1, i++) {
2055	if (j & 1) {
2056	ieps++;
2057	dval(&u) *= bigtens[i];
2058	}
2059	}
2060	}
2061	if (k_check && dval(&u) < 1. && ilim > 0) {
2062	if (ilim1 <= 0)
2063	goto fast_failed;
2064	ilim = ilim1;
2065	k--;
2066	dval(&u) *= 10.;
2067	ieps++;
2068	}
2069	dval(&eps) = (ieps * dval(&u)) + 7.;
2070	word0(&eps) -= (P - 1) * Exp_msk1;
2071	if (ilim == 0) {
2072	S.clear();
2073	mhi.clear();
2074	dval(&u) -= 5.;
2075	if (dval(&u) > dval(&eps))
2076	goto one_digit;
2077	if (dval(&u) < -dval(&eps))
2078	goto no_digits;
2079	goto fast_failed;
2080	}
2081	#ifndef No_leftright
2082	if (leftright) {
2083	/* Use Steele & White method of only
2084	* generating digits needed.
2085	*/
2086	dval(&eps) = (0.5 / tens[ilim - 1]) - dval(&eps);
2087	for (i = 0;;) {
2088	L = (long int)dval(&u);
2089	dval(&u) -= L;
2090	*s++ = '0' + (int)L;
2091	if (dval(&u) < dval(&eps))
2092	goto ret;
2093	if (1. - dval(&u) < dval(&eps))
2094	goto bump_up;
2095	if (++i >= ilim)
2096	break;
2097	dval(&eps) *= 10.;
2098	dval(&u) *= 10.;
2099	}
2100	} else {
2101	#endif
2102	/* Generate ilim digits, then fix them up. */
2103	dval(&eps) *= tens[ilim - 1];
2104	for (i = 1;; i++, dval(&u) *= 10.) {
2105	L = (int32_t)(dval(&u));
2106	if (!(dval(&u) -= L))
2107	ilim = i;
2108	*s++ = '0' + (int)L;
2109	if (i == ilim) {
2110	if (dval(&u) > 0.5 + dval(&eps))
2111	goto bump_up;
2112	else if (dval(&u) < 0.5 - dval(&eps)) {
2113	while (*--s == '0') { }
2114	s++;
2115	goto ret;
2116	}
2117	break;
2118	}
2119	}
2120	#ifndef No_leftright
2121	}
2122	#endif
2123	fast_failed:
2124	s = s0;
2125	dval(&u) = dval(&d2);
2126	k = k0;
2127	ilim = ilim0;
2128	}
2129
2130	/* Do we have a "small" integer? */
2131
2132	if (be >= 0 && k <= Int_max) {
2133	/* Yes. */
2134	ds = tens[k];
2135	if (ndigits < 0 && ilim <= 0) {
2136	S.clear();
2137	mhi.clear();
2138	if (ilim < 0 || dval(&u) <= 5 * ds)
2139	goto no_digits;
2140	goto one_digit;
2141	}
2142	for (i = 1;; i++, dval(&u) *= 10.) {
2143	L = (int32_t)(dval(&u) / ds);
2144	dval(&u) -= L * ds;
2145	#ifdef Check_FLT_ROUNDS
2146	/* If FLT_ROUNDS == 2, L will usually be high by 1 */
2147	if (dval(&u) < 0) {
2148	L--;
2149	dval(&u) += ds;
2150	}
2151	#endif
2152	*s++ = '0' + (int)L;
2153	if (!dval(&u)) {
2154	#ifdef SET_INEXACT
2155	inexact = 0;
2156	#endif
2157	break;
2158	}
2159	if (i == ilim) {
2160	dval(&u) += dval(&u);
2161	if (dval(&u) > ds || (dval(&u) == ds && (L & 1))) {
2162	bump_up:
2163	while (*--s == '9')
2164	if (s == s0) {
2165	k++;
2166	*s = '0';
2167	break;
2168	}
2169	++*s++;
2170	}
2171	break;
2172	}
2173	}
2174	goto ret;
2175	}
2176
2177	m2 = b2;
2178	m5 = b5;
2179	mhi.clear();
2180	mlo.clear();
2181	if (leftright) {
2182	i =
2183	#ifndef Sudden_Underflow
2184	denorm ? be + (Bias + (P - 1) - 1 + 1) :
2185	#endif
2186	1 + P - bbits;
2187	b2 += i;
2188	s2 += i;
2189	i2b(b&: mhi, i: 1);
2190	}
2191	if (m2 > 0 && s2 > 0) {
2192	i = m2 < s2 ? m2 : s2;
2193	b2 -= i;
2194	m2 -= i;
2195	s2 -= i;
2196	}
2197	if (b5 > 0) {
2198	if (leftright) {
2199	if (m5 > 0) {
2200	pow5mult(b&: mhi, k: m5);
2201	mult(aRef&: b, bRef: mhi);
2202	}
2203	if ((j = b5 - m5))
2204	pow5mult(b, k: j);
2205	} else
2206	pow5mult(b, k: b5);
2207	}
2208	i2b(b&: S, i: 1);
2209	if (s5 > 0)
2210	pow5mult(b&: S, k: s5);
2211
2212	/* Check for special case that d is a normalized power of 2. */
2213
2214	spec_case = 0;
2215	if (!word1(&u) && !(word0(&u) & Bndry_mask)
2216	#ifndef Sudden_Underflow
2217	&& word0(&u) & (Exp_mask & ~Exp_msk1)
2218	#endif
2219	) {
2220	/* The special case */
2221	b2 += Log2P;
2222	s2 += Log2P;
2223	spec_case = 1;
2224	}
2225
2226	/* Arrange for convenient computation of quotients:
2227	* shift left if necessary so divisor has 4 leading 0 bits.
2228	*
2229	* Perhaps we should just compute leading 28 bits of S once
2230	* and for all and pass them and a shift to quorem, so it
2231	* can do shifts and ors to compute the numerator for q.
2232	*/
2233	#ifdef Pack_32
2234	if ((i = ((s5 ? 32 - hi0bits(x: S.words()[S.size() - 1]) : 1) + s2) & 0x1f))
2235	i = 32 - i;
2236	#else
2237	if ((i = ((s5 ? 32 - hi0bits(S.words()[S.size() - 1]) : 1) + s2) & 0xf))
2238	i = 16 - i;
2239	#endif
2240	if (i > 4) {
2241	i -= 4;
2242	b2 += i;
2243	m2 += i;
2244	s2 += i;
2245	} else if (i < 4) {
2246	i += 28;
2247	b2 += i;
2248	m2 += i;
2249	s2 += i;
2250	}
2251	if (b2 > 0)
2252	lshift(b, k: b2);
2253	if (s2 > 0)
2254	lshift(b&: S, k: s2);
2255	if (k_check) {
2256	if (cmp(a: b,b: S) < 0) {
2257	k--;
2258	multadd(b, m: 10, a: 0); /* we botched the k estimate */
2259	if (leftright)
2260	multadd(b&: mhi, m: 10, a: 0);
2261	ilim = ilim1;
2262	}
2263	}
2264
2265	if (leftright) {
2266	if (m2 > 0)
2267	lshift(b&: mhi, k: m2);
2268
2269	/* Compute mlo -- check for special case
2270	* that d is a normalized power of 2.
2271	*/
2272
2273	mlo = mhi;
2274	if (spec_case) {
2275	mhi = mlo;
2276	lshift(b&: mhi, Log2P);
2277	}
2278
2279	for (i = 1;;i++) {
2280	dig = quorem(b,S) + '0';
2281	/* Do we yet have the shortest decimal string
2282	* that will round to d?
2283	*/
2284	j = cmp(a: b, b: mlo);
2285	diff(c&: delta, aRef: S, bRef: mhi);
2286	j1 = delta.sign ? 1 : cmp(a: b, b: delta);
2287	if (j1 == 0 && !(word1(&u) & 1)) {
2288	if (dig == '9')
2289	goto round_9_up;
2290	if (j > 0)
2291	dig++;
2292	#ifdef SET_INEXACT
2293	else if (!b->x[0] && b->wds <= 1)
2294	inexact = 0;
2295	#endif
2296	*s++ = dig;
2297	goto ret;
2298	}
2299	if (j < 0 || (j == 0 && !(word1(&u) & 1))) {
2300	if (!b.words()[0] && b.size() <= 1) {
2301	#ifdef SET_INEXACT
2302	inexact = 0;
2303	#endif
2304	goto accept_dig;
2305	}
2306	if (j1 > 0) {
2307	lshift(b, k: 1);
2308	j1 = cmp(a: b, b: S);
2309	if ((j1 > 0 || (j1 == 0 && (dig & 1))) && dig++ == '9')
2310	goto round_9_up;
2311	}
2312	accept_dig:
2313	*s++ = dig;
2314	goto ret;
2315	}
2316	if (j1 > 0) {
2317	if (dig == '9') { /* possible if i == 1 */
2318	round_9_up:
2319	*s++ = '9';
2320	goto roundoff;
2321	}
2322	*s++ = dig + 1;
2323	goto ret;
2324	}
2325	*s++ = dig;
2326	if (i == ilim)
2327	break;
2328	multadd(b, m: 10, a: 0);
2329	multadd(b&: mlo, m: 10, a: 0);
2330	multadd(b&: mhi, m: 10, a: 0);
2331	}
2332	} else
2333	for (i = 1;; i++) {
2334	*s++ = dig = quorem(b,S) + '0';
2335	if (!b.words()[0] && b.size() <= 1) {
2336	#ifdef SET_INEXACT
2337	inexact = 0;
2338	#endif
2339	goto ret;
2340	}
2341	if (i >= ilim)
2342	break;
2343	multadd(b, m: 10, a: 0);
2344	}
2345
2346	/* Round off last digit */
2347
2348	lshift(b, k: 1);
2349	j = cmp(a: b, b: S);
2350	if (j > 0 || (j == 0 && (dig & 1))) {
2351	roundoff:
2352	while (*--s == '9')
2353	if (s == s0) {
2354	k++;
2355	*s++ = '1';
2356	goto ret;
2357	}
2358	++*s++;
2359	} else {
2360	while (*--s == '0') { }
2361	s++;
2362	}
2363	goto ret;
2364	no_digits:
2365	k = -1 - ndigits;
2366	goto ret;
2367	one_digit:
2368	*s++ = '1';
2369	k++;
2370	goto ret;
2371	ret:
2372	#ifdef SET_INEXACT
2373	if (inexact) {
2374	if (!oldinexact) {
2375	word0(&u) = Exp_1 + (70 << Exp_shift);
2376	word1(&u) = 0;
2377	dval(&u) += 1.;
2378	}
2379	} else if (!oldinexact)
2380	clear_inexact();
2381	#endif
2382	*s = 0;
2383	*decpt = k + 1;
2384	if (rve)
2385	*rve = s;
2386	}
2387
2388	static ALWAYS_INLINE void append(char*& next, const char* src, unsigned size)
2389	{
2390	for (unsigned i = 0; i < size; ++i)
2391	*next++ = *src++;
2392	}
2393
2394	void doubleToStringInJavaScriptFormat(double d, DtoaBuffer buffer, unsigned* resultLength)
2395	{
2396	ASSERT(buffer);
2397
2398	// avoid ever printing -NaN, in JS conceptually there is only one NaN value
2399	if (std::isnan(x: d)) {
2400	append(next&: buffer, src: "NaN", size: 3);
2401	if (resultLength)
2402	*resultLength = 3;
2403	return;
2404	}
2405	// -0 -> "0"
2406	if (!d) {
2407	buffer[0] = '0';
2408	if (resultLength)
2409	*resultLength = 1;
2410	return;
2411	}
2412
2413	int decimalPoint;
2414	int sign;
2415
2416	DtoaBuffer result;
2417	char* resultEnd = 0;
2418	WTF::dtoa(result, dd: d, ndigits: 0, decpt: &decimalPoint, sign: &sign, rve: &resultEnd);
2419	int length = resultEnd - result;
2420
2421	char* next = buffer;
2422	if (sign)
2423	*next++ = '-';
2424
2425	if (decimalPoint <= 0 && decimalPoint > -6) {
2426	*next++ = '0';
2427	*next++ = '.';
2428	for (int j = decimalPoint; j < 0; j++)
2429	*next++ = '0';
2430	append(next, src: result, size: length);
2431	} else if (decimalPoint <= 21 && decimalPoint > 0) {
2432	if (length <= decimalPoint) {
2433	append(next, src: result, size: length);
2434	for (int j = 0; j < decimalPoint - length; j++)
2435	*next++ = '0';
2436	} else {
2437	append(next, src: result, size: decimalPoint);
2438	*next++ = '.';
2439	append(next, src: result + decimalPoint, size: length - decimalPoint);
2440	}
2441	} else if (result[0] < '0' || result[0] > '9')
2442	append(next, src: result, size: length);
2443	else {
2444	*next++ = result[0];
2445	if (length > 1) {
2446	*next++ = '.';
2447	append(next, src: result + 1, size: length - 1);
2448	}
2449
2450	*next++ = 'e';
2451	*next++ = (decimalPoint >= 0) ? '+' : '-';
2452	// decimalPoint can't be more than 3 digits decimal given the
2453	// nature of float representation
2454	int exponential = decimalPoint - 1;
2455	if (exponential < 0)
2456	exponential = -exponential;
2457	if (exponential >= 100)
2458	*next++ = static_cast<char>('0' + exponential / 100);
2459	if (exponential >= 10)
2460	*next++ = static_cast<char>('0' + (exponential % 100) / 10);
2461	*next++ = static_cast<char>('0' + exponential % 10);
2462	}
2463	if (resultLength)
2464	*resultLength = next - buffer;
2465	}
2466
2467	} // namespace WTF
2468