• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1 /****************************************************************
2  *
3  * The author of this software is David M. Gay.
4  *
5  * Copyright (c) 1991, 2000, 2001 by Lucent Technologies.
6  *
7  * Permission to use, copy, modify, and distribute this software for any
8  * purpose without fee is hereby granted, provided that this entire notice
9  * is included in all copies of any software which is or includes a copy
10  * or modification of this software and in all copies of the supporting
11  * documentation for such software.
12  *
13  * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
14  * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY
15  * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
16  * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
17  *
18  ***************************************************************/
19 
20 /* Please send bug reports to David M. Gay (dmg at acm dot org,
21  * with " at " changed at "@" and " dot " changed to ".").	*/
22 
23 /* On a machine with IEEE extended-precision registers, it is
24  * necessary to specify double-precision (53-bit) rounding precision
25  * before invoking strtod or dtoa.  If the machine uses (the equivalent
26  * of) Intel 80x87 arithmetic, the call
27  *	_control87(PC_53, MCW_PC);
28  * does this with many compilers.  Whether this or another call is
29  * appropriate depends on the compiler; for this to work, it may be
30  * necessary to #include "float.h" or another system-dependent header
31  * file.
32  */
33 
34 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
35  *
36  * This strtod returns a nearest machine number to the input decimal
37  * string (or sets errno to ERANGE).  With IEEE arithmetic, ties are
38  * broken by the IEEE round-even rule.  Otherwise ties are broken by
39  * biased rounding (add half and chop).
40  *
41  * Inspired loosely by William D. Clinger's paper "How to Read Floating
42  * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
43  *
44  * Modifications:
45  *
46  *	1. We only require IEEE, IBM, or VAX double-precision
47  *		arithmetic (not IEEE double-extended).
48  *	2. We get by with floating-point arithmetic in a case that
49  *		Clinger missed -- when we're computing d * 10^n
50  *		for a small integer d and the integer n is not too
51  *		much larger than 22 (the maximum integer k for which
52  *		we can represent 10^k exactly), we may be able to
53  *		compute (d*10^k) * 10^(e-k) with just one roundoff.
54  *	3. Rather than a bit-at-a-time adjustment of the binary
55  *		result in the hard case, we use floating-point
56  *		arithmetic to determine the adjustment to within
57  *		one bit; only in really hard cases do we need to
58  *		compute a second residual.
59  *	4. Because of 3., we don't need a large table of powers of 10
60  *		for ten-to-e (just some small tables, e.g. of 10^k
61  *		for 0 <= k <= 22).
62  */
63 
64 /*
65  * #define IEEE_8087 for IEEE-arithmetic machines where the least
66  *	significant byte has the lowest address.
67  * #define IEEE_MC68k for IEEE-arithmetic machines where the most
68  *	significant byte has the lowest address.
69  * #define Long int on machines with 32-bit ints and 64-bit longs.
70  * #define IBM for IBM mainframe-style floating-point arithmetic.
71  * #define VAX for VAX-style floating-point arithmetic (D_floating).
72  * #define No_leftright to omit left-right logic in fast floating-point
73  *	computation of dtoa.
74  * #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
75  *	and strtod and dtoa should round accordingly.  Unless Trust_FLT_ROUNDS
76  *	is also #defined, fegetround() will be queried for the rounding mode.
77  *	Note that both FLT_ROUNDS and fegetround() are specified by the C99
78  *	standard (and are specified to be consistent, with fesetround()
79  *	affecting the value of FLT_ROUNDS), but that some (Linux) systems
80  *	do not work correctly in this regard, so using fegetround() is more
81  *	portable than using FLT_FOUNDS directly.
82  * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
83  *	and Honor_FLT_ROUNDS is not #defined.
84  * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
85  *	that use extended-precision instructions to compute rounded
86  *	products and quotients) with IBM.
87  * #define ROUND_BIASED for IEEE-format with biased rounding.
88  * #define Inaccurate_Divide for IEEE-format with correctly rounded
89  *	products but inaccurate quotients, e.g., for Intel i860.
90  * #define NO_LONG_LONG on machines that do not have a "long long"
91  *	integer type (of >= 64 bits).  On such machines, you can
92  *	#define Just_16 to store 16 bits per 32-bit Long when doing
93  *	high-precision integer arithmetic.  Whether this speeds things
94  *	up or slows things down depends on the machine and the number
95  *	being converted.  If long long is available and the name is
96  *	something other than "long long", #define Llong to be the name,
97  *	and if "unsigned Llong" does not work as an unsigned version of
98  *	Llong, #define #ULLong to be the corresponding unsigned type.
99  * #define KR_headers for old-style C function headers.
100  * #define Bad_float_h if your system lacks a float.h or if it does not
101  *	define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
102  *	FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
103  * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
104  *	if memory is available and otherwise does something you deem
105  *	appropriate.  If MALLOC is undefined, malloc will be invoked
106  *	directly -- and assumed always to succeed.  Similarly, if you
107  *	want something other than the system's free() to be called to
108  *	recycle memory acquired from MALLOC, #define FREE to be the
109  *	name of the alternate routine.  (FREE or free is only called in
110  *	pathological cases, e.g., in a dtoa call after a dtoa return in
111  *	mode 3 with thousands of digits requested.)
112  * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making
113  *	memory allocations from a private pool of memory when possible.
114  *	When used, the private pool is PRIVATE_MEM bytes long:  2304 bytes,
115  *	unless #defined to be a different length.  This default length
116  *	suffices to get rid of MALLOC calls except for unusual cases,
117  *	such as decimal-to-binary conversion of a very long string of
118  *	digits.  The longest string dtoa can return is about 751 bytes
119  *	long.  For conversions by strtod of strings of 800 digits and
120  *	all dtoa conversions in single-threaded executions with 8-byte
121  *	pointers, PRIVATE_MEM >= 7400 appears to suffice; with 4-byte
122  *	pointers, PRIVATE_MEM >= 7112 appears adequate.
123  * #define NO_INFNAN_CHECK if you do not wish to have INFNAN_CHECK
124  *	#defined automatically on IEEE systems.  On such systems,
125  *	when INFNAN_CHECK is #defined, strtod checks
126  *	for Infinity and NaN (case insensitively).  On some systems
127  *	(e.g., some HP systems), it may be necessary to #define NAN_WORD0
128  *	appropriately -- to the most significant word of a quiet NaN.
129  *	(On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.)
130  *	When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined,
131  *	strtod also accepts (case insensitively) strings of the form
132  *	NaN(x), where x is a string of hexadecimal digits and spaces;
133  *	if there is only one string of hexadecimal digits, it is taken
134  *	for the 52 fraction bits of the resulting NaN; if there are two
135  *	or more strings of hex digits, the first is for the high 20 bits,
136  *	the second and subsequent for the low 32 bits, with intervening
137  *	white space ignored; but if this results in none of the 52
138  *	fraction bits being on (an IEEE Infinity symbol), then NAN_WORD0
139  *	and NAN_WORD1 are used instead.
140  * #define MULTIPLE_THREADS if the system offers preemptively scheduled
141  *	multiple threads.  In this case, you must provide (or suitably
142  *	#define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed
143  *	by FREE_DTOA_LOCK(n) for n = 0 or 1.  (The second lock, accessed
144  *	in pow5mult, ensures lazy evaluation of only one copy of high
145  *	powers of 5; omitting this lock would introduce a small
146  *	probability of wasting memory, but would otherwise be harmless.)
147  *	You must also invoke freedtoa(s) to free the value s returned by
148  *	dtoa.  You may do so whether or not MULTIPLE_THREADS is #defined.
149  * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that
150  *	avoids underflows on inputs whose result does not underflow.
151  *	If you #define NO_IEEE_Scale on a machine that uses IEEE-format
152  *	floating-point numbers and flushes underflows to zero rather
153  *	than implementing gradual underflow, then you must also #define
154  *	Sudden_Underflow.
155  * #define USE_LOCALE to use the current locale's decimal_point value.
156  * #define SET_INEXACT if IEEE arithmetic is being used and extra
157  *	computation should be done to set the inexact flag when the
158  *	result is inexact and avoid setting inexact when the result
159  *	is exact.  In this case, dtoa.c must be compiled in
160  *	an environment, perhaps provided by #include "dtoa.c" in a
161  *	suitable wrapper, that defines two functions,
162  *		int get_inexact(void);
163  *		void clear_inexact(void);
164  *	such that get_inexact() returns a nonzero value if the
165  *	inexact bit is already set, and clear_inexact() sets the
166  *	inexact bit to 0.  When SET_INEXACT is #defined, strtod
167  *	also does extra computations to set the underflow and overflow
168  *	flags when appropriate (i.e., when the result is tiny and
169  *	inexact or when it is a numeric value rounded to +-infinity).
170  * #define NO_ERRNO if strtod should not assign errno = ERANGE when
171  *	the result overflows to +-Infinity or underflows to 0.
172  * #define NO_HEX_FP to omit recognition of hexadecimal floating-point
173  *	values by strtod.
174  * #define NO_STRTOD_BIGCOMP (on IEEE-arithmetic systems only for now)
175  *	to disable logic for "fast" testing of very long input strings
176  *	to strtod.  This testing proceeds by initially truncating the
177  *	input string, then if necessary comparing the whole string with
178  *	a decimal expansion to decide close cases. This logic is only
179  *	used for input more than STRTOD_DIGLIM digits long (default 40).
180  */
181 
182 #define IEEE_8087
183 #define NO_HEX_FP
184 
185 #ifndef Long
186 #if __LP64__
187 #define Long int
188 #else
189 #define Long long
190 #endif
191 #endif
192 #ifndef ULong
193 typedef unsigned Long ULong;
194 #endif
195 
196 #ifdef DEBUG
197 #include "stdio.h"
198 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
199 #endif
200 
201 #include "stdlib.h"
202 #include "string.h"
203 
204 #ifdef USE_LOCALE
205 #include "locale.h"
206 #endif
207 
208 #ifdef Honor_FLT_ROUNDS
209 #ifndef Trust_FLT_ROUNDS
210 #include <fenv.h>
211 #endif
212 #endif
213 
214 #ifdef MALLOC
215 #ifdef KR_headers
216 extern char *MALLOC();
217 #else
218 extern void *MALLOC(size_t);
219 #endif
220 #else
221 #define MALLOC malloc
222 #endif
223 
224 #ifndef Omit_Private_Memory
225 #ifndef PRIVATE_MEM
226 #define PRIVATE_MEM 2304
227 #endif
228 #define PRIVATE_mem ((unsigned)((PRIVATE_MEM+sizeof(double)-1)/sizeof(double)))
229 static double private_mem[PRIVATE_mem], *pmem_next = private_mem;
230 #endif
231 
232 #undef IEEE_Arith
233 #undef Avoid_Underflow
234 #ifdef IEEE_MC68k
235 #define IEEE_Arith
236 #endif
237 #ifdef IEEE_8087
238 #define IEEE_Arith
239 #endif
240 
241 #ifdef IEEE_Arith
242 #ifndef NO_INFNAN_CHECK
243 #undef INFNAN_CHECK
244 #define INFNAN_CHECK
245 #endif
246 #else
247 #undef INFNAN_CHECK
248 #define NO_STRTOD_BIGCOMP
249 #endif
250 
251 #include "errno.h"
252 
253 #ifdef Bad_float_h
254 
255 #ifdef IEEE_Arith
256 #define DBL_DIG 15
257 #define DBL_MAX_10_EXP 308
258 #define DBL_MAX_EXP 1024
259 #define FLT_RADIX 2
260 #endif /*IEEE_Arith*/
261 
262 #ifdef IBM
263 #define DBL_DIG 16
264 #define DBL_MAX_10_EXP 75
265 #define DBL_MAX_EXP 63
266 #define FLT_RADIX 16
267 #define DBL_MAX 7.2370055773322621e+75
268 #endif
269 
270 #ifdef VAX
271 #define DBL_DIG 16
272 #define DBL_MAX_10_EXP 38
273 #define DBL_MAX_EXP 127
274 #define FLT_RADIX 2
275 #define DBL_MAX 1.7014118346046923e+38
276 #endif
277 
278 #ifndef LONG_MAX
279 #define LONG_MAX 2147483647
280 #endif
281 
282 #else /* ifndef Bad_float_h */
283 #include "float.h"
284 #endif /* Bad_float_h */
285 
286 #ifndef __MATH_H__
287 #include "math.h"
288 #endif
289 
290 namespace dmg_fp {
291 
292 #ifndef CONST
293 #ifdef KR_headers
294 #define CONST /* blank */
295 #else
296 #define CONST const
297 #endif
298 #endif
299 
300 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
301 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
302 #endif
303 
304 typedef union { double d; ULong L[2]; } U;
305 
306 #ifdef IEEE_8087
307 #define word0(x) (x)->L[1]
308 #define word1(x) (x)->L[0]
309 #else
310 #define word0(x) (x)->L[0]
311 #define word1(x) (x)->L[1]
312 #endif
313 #define dval(x) (x)->d
314 
315 #ifndef STRTOD_DIGLIM
316 #define STRTOD_DIGLIM 40
317 #endif
318 
319 #ifdef DIGLIM_DEBUG
320 extern int strtod_diglim;
321 #else
322 #define strtod_diglim STRTOD_DIGLIM
323 #endif
324 
325 /* The following definition of Storeinc is appropriate for MIPS processors.
326  * An alternative that might be better on some machines is
327  * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
328  */
329 #if defined(IEEE_8087) + defined(VAX)
330 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
331 ((unsigned short *)a)[0] = (unsigned short)c, a++)
332 #else
333 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
334 ((unsigned short *)a)[1] = (unsigned short)c, a++)
335 #endif
336 
337 /* #define P DBL_MANT_DIG */
338 /* Ten_pmax = floor(P*log(2)/log(5)) */
339 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
340 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
341 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
342 
343 #ifdef IEEE_Arith
344 #define Exp_shift  20
345 #define Exp_shift1 20
346 #define Exp_msk1    0x100000
347 #define Exp_msk11   0x100000
348 #define Exp_mask  0x7ff00000
349 #define P 53
350 #define Nbits 53
351 #define Bias 1023
352 #define Emax 1023
353 #define Emin (-1022)
354 #define Exp_1  0x3ff00000
355 #define Exp_11 0x3ff00000
356 #define Ebits 11
357 #define Frac_mask  0xfffff
358 #define Frac_mask1 0xfffff
359 #define Ten_pmax 22
360 #define Bletch 0x10
361 #define Bndry_mask  0xfffff
362 #define Bndry_mask1 0xfffff
363 #define LSB 1
364 #define Sign_bit 0x80000000
365 #define Log2P 1
366 #define Tiny0 0
367 #define Tiny1 1
368 #define Quick_max 14
369 #define Int_max 14
370 #ifndef NO_IEEE_Scale
371 #define Avoid_Underflow
372 #ifdef Flush_Denorm	/* debugging option */
373 #undef Sudden_Underflow
374 #endif
375 #endif
376 
377 #ifndef Flt_Rounds
378 #ifdef FLT_ROUNDS
379 #define Flt_Rounds FLT_ROUNDS
380 #else
381 #define Flt_Rounds 1
382 #endif
383 #endif /*Flt_Rounds*/
384 
385 #ifdef Honor_FLT_ROUNDS
386 #undef Check_FLT_ROUNDS
387 #define Check_FLT_ROUNDS
388 #else
389 #define Rounding Flt_Rounds
390 #endif
391 
392 #else /* ifndef IEEE_Arith */
393 #undef Check_FLT_ROUNDS
394 #undef Honor_FLT_ROUNDS
395 #undef SET_INEXACT
396 #undef  Sudden_Underflow
397 #define Sudden_Underflow
398 #ifdef IBM
399 #undef Flt_Rounds
400 #define Flt_Rounds 0
401 #define Exp_shift  24
402 #define Exp_shift1 24
403 #define Exp_msk1   0x1000000
404 #define Exp_msk11  0x1000000
405 #define Exp_mask  0x7f000000
406 #define P 14
407 #define Nbits 56
408 #define Bias 65
409 #define Emax 248
410 #define Emin (-260)
411 #define Exp_1  0x41000000
412 #define Exp_11 0x41000000
413 #define Ebits 8	/* exponent has 7 bits, but 8 is the right value in b2d */
414 #define Frac_mask  0xffffff
415 #define Frac_mask1 0xffffff
416 #define Bletch 4
417 #define Ten_pmax 22
418 #define Bndry_mask  0xefffff
419 #define Bndry_mask1 0xffffff
420 #define LSB 1
421 #define Sign_bit 0x80000000
422 #define Log2P 4
423 #define Tiny0 0x100000
424 #define Tiny1 0
425 #define Quick_max 14
426 #define Int_max 15
427 #else /* VAX */
428 #undef Flt_Rounds
429 #define Flt_Rounds 1
430 #define Exp_shift  23
431 #define Exp_shift1 7
432 #define Exp_msk1    0x80
433 #define Exp_msk11   0x800000
434 #define Exp_mask  0x7f80
435 #define P 56
436 #define Nbits 56
437 #define Bias 129
438 #define Emax 126
439 #define Emin (-129)
440 #define Exp_1  0x40800000
441 #define Exp_11 0x4080
442 #define Ebits 8
443 #define Frac_mask  0x7fffff
444 #define Frac_mask1 0xffff007f
445 #define Ten_pmax 24
446 #define Bletch 2
447 #define Bndry_mask  0xffff007f
448 #define Bndry_mask1 0xffff007f
449 #define LSB 0x10000
450 #define Sign_bit 0x8000
451 #define Log2P 1
452 #define Tiny0 0x80
453 #define Tiny1 0
454 #define Quick_max 15
455 #define Int_max 15
456 #endif /* IBM, VAX */
457 #endif /* IEEE_Arith */
458 
459 #ifndef IEEE_Arith
460 #define ROUND_BIASED
461 #endif
462 
463 #ifdef RND_PRODQUOT
464 #define rounded_product(a,b) a = rnd_prod(a, b)
465 #define rounded_quotient(a,b) a = rnd_quot(a, b)
466 #ifdef KR_headers
467 extern double rnd_prod(), rnd_quot();
468 #else
469 extern double rnd_prod(double, double), rnd_quot(double, double);
470 #endif
471 #else
472 #define rounded_product(a,b) a *= b
473 #define rounded_quotient(a,b) a /= b
474 #endif
475 
476 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
477 #define Big1 0xffffffff
478 
479 #ifndef Pack_32
480 #define Pack_32
481 #endif
482 
483 typedef struct BCinfo BCinfo;
484  struct
485 BCinfo { int dp0, dp1, dplen, dsign, e0, inexact, nd, nd0, rounding, scale, uflchk; };
486 
487 #ifdef KR_headers
488 #define FFFFFFFF ((((unsigned long)0xffff)<<16)|(unsigned long)0xffff)
489 #else
490 #define FFFFFFFF 0xffffffffUL
491 #endif
492 
493 #ifdef NO_LONG_LONG
494 #undef ULLong
495 #ifdef Just_16
496 #undef Pack_32
497 /* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
498  * This makes some inner loops simpler and sometimes saves work
499  * during multiplications, but it often seems to make things slightly
500  * slower.  Hence the default is now to store 32 bits per Long.
501  */
502 #endif
503 #else	/* long long available */
504 #ifndef Llong
505 #define Llong long long
506 #endif
507 #ifndef ULLong
508 #define ULLong unsigned Llong
509 #endif
510 #endif /* NO_LONG_LONG */
511 
512 #ifndef MULTIPLE_THREADS
513 #define ACQUIRE_DTOA_LOCK(n)	/*nothing*/
514 #define FREE_DTOA_LOCK(n)	/*nothing*/
515 #endif
516 
517 #define Kmax 7
518 
519 double strtod(const char *s00, char **se);
520 char *dtoa(double d, int mode, int ndigits,
521 			int *decpt, int *sign, char **rve);
522 
523  struct
524 Bigint {
525 	struct Bigint *next;
526 	int k, maxwds, sign, wds;
527 	ULong x[1];
528 	};
529 
530  typedef struct Bigint Bigint;
531 
532  static Bigint *freelist[Kmax+1];
533 
534  static Bigint *
535 Balloc
536 #ifdef KR_headers
537 	(k) int k;
538 #else
539 	(int k)
540 #endif
541 {
542 	int x;
543 	Bigint *rv;
544 #ifndef Omit_Private_Memory
545 	unsigned int len;
546 #endif
547 
548 	ACQUIRE_DTOA_LOCK(0);
549 	/* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */
550 	/* but this case seems very unlikely. */
551 	if (k <= Kmax && (rv = freelist[k]))
552 		freelist[k] = rv->next;
553 	else {
554 		x = 1 << k;
555 #ifdef Omit_Private_Memory
556 		rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong));
557 #else
558 		len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
559 			/sizeof(double);
560 		if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) {
561 			rv = (Bigint*)pmem_next;
562 			pmem_next += len;
563 			}
564 		else
565 			rv = (Bigint*)MALLOC(len*sizeof(double));
566 #endif
567 		rv->k = k;
568 		rv->maxwds = x;
569 		}
570 	FREE_DTOA_LOCK(0);
571 	rv->sign = rv->wds = 0;
572 	return rv;
573 	}
574 
575  static void
576 Bfree
577 #ifdef KR_headers
578 	(v) Bigint *v;
579 #else
580 	(Bigint *v)
581 #endif
582 {
583 	if (v) {
584 		if (v->k > Kmax)
585 #ifdef FREE
586 			FREE((void*)v);
587 #else
588 			free((void*)v);
589 #endif
590 		else {
591 			ACQUIRE_DTOA_LOCK(0);
592 			v->next = freelist[v->k];
593 			freelist[v->k] = v;
594 			FREE_DTOA_LOCK(0);
595 			}
596 		}
597 	}
598 
599 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
600 y->wds*sizeof(Long) + 2*sizeof(int))
601 
602  static Bigint *
603 multadd
604 #ifdef KR_headers
605 	(b, m, a) Bigint *b; int m, a;
606 #else
607 	(Bigint *b, int m, int a)	/* multiply by m and add a */
608 #endif
609 {
610 	int i, wds;
611 #ifdef ULLong
612 	ULong *x;
613 	ULLong carry, y;
614 #else
615 	ULong carry, *x, y;
616 #ifdef Pack_32
617 	ULong xi, z;
618 #endif
619 #endif
620 	Bigint *b1;
621 
622 	wds = b->wds;
623 	x = b->x;
624 	i = 0;
625 	carry = a;
626 	do {
627 #ifdef ULLong
628 		y = *x * (ULLong)m + carry;
629 		carry = y >> 32;
630 		*x++ = y & FFFFFFFF;
631 #else
632 #ifdef Pack_32
633 		xi = *x;
634 		y = (xi & 0xffff) * m + carry;
635 		z = (xi >> 16) * m + (y >> 16);
636 		carry = z >> 16;
637 		*x++ = (z << 16) + (y & 0xffff);
638 #else
639 		y = *x * m + carry;
640 		carry = y >> 16;
641 		*x++ = y & 0xffff;
642 #endif
643 #endif
644 		}
645 		while(++i < wds);
646 	if (carry) {
647 		if (wds >= b->maxwds) {
648 			b1 = Balloc(b->k+1);
649 			Bcopy(b1, b);
650 			Bfree(b);
651 			b = b1;
652 			}
653 		b->x[wds++] = carry;
654 		b->wds = wds;
655 		}
656 	return b;
657 	}
658 
659  static Bigint *
660 s2b
661 #ifdef KR_headers
662 	(s, nd0, nd, y9, dplen) CONST char *s; int nd0, nd, dplen; ULong y9;
663 #else
664 	(CONST char *s, int nd0, int nd, ULong y9, int dplen)
665 #endif
666 {
667 	Bigint *b;
668 	int i, k;
669 	Long x, y;
670 
671 	x = (nd + 8) / 9;
672 	for(k = 0, y = 1; x > y; y <<= 1, k++) ;
673 #ifdef Pack_32
674 	b = Balloc(k);
675 	b->x[0] = y9;
676 	b->wds = 1;
677 #else
678 	b = Balloc(k+1);
679 	b->x[0] = y9 & 0xffff;
680 	b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
681 #endif
682 
683 	i = 9;
684 	if (9 < nd0) {
685 		s += 9;
686 		do b = multadd(b, 10, *s++ - '0');
687 			while(++i < nd0);
688 		s += dplen;
689 		}
690 	else
691 		s += dplen + 9;
692 	for(; i < nd; i++)
693 		b = multadd(b, 10, *s++ - '0');
694 	return b;
695 	}
696 
697  static int
698 hi0bits
699 #ifdef KR_headers
700 	(x) ULong x;
701 #else
702 	(ULong x)
703 #endif
704 {
705 	int k = 0;
706 
707 	if (!(x & 0xffff0000)) {
708 		k = 16;
709 		x <<= 16;
710 		}
711 	if (!(x & 0xff000000)) {
712 		k += 8;
713 		x <<= 8;
714 		}
715 	if (!(x & 0xf0000000)) {
716 		k += 4;
717 		x <<= 4;
718 		}
719 	if (!(x & 0xc0000000)) {
720 		k += 2;
721 		x <<= 2;
722 		}
723 	if (!(x & 0x80000000)) {
724 		k++;
725 		if (!(x & 0x40000000))
726 			return 32;
727 		}
728 	return k;
729 	}
730 
731  static int
732 lo0bits
733 #ifdef KR_headers
734 	(y) ULong *y;
735 #else
736 	(ULong *y)
737 #endif
738 {
739 	int k;
740 	ULong x = *y;
741 
742 	if (x & 7) {
743 		if (x & 1)
744 			return 0;
745 		if (x & 2) {
746 			*y = x >> 1;
747 			return 1;
748 			}
749 		*y = x >> 2;
750 		return 2;
751 		}
752 	k = 0;
753 	if (!(x & 0xffff)) {
754 		k = 16;
755 		x >>= 16;
756 		}
757 	if (!(x & 0xff)) {
758 		k += 8;
759 		x >>= 8;
760 		}
761 	if (!(x & 0xf)) {
762 		k += 4;
763 		x >>= 4;
764 		}
765 	if (!(x & 0x3)) {
766 		k += 2;
767 		x >>= 2;
768 		}
769 	if (!(x & 1)) {
770 		k++;
771 		x >>= 1;
772 		if (!x)
773 			return 32;
774 		}
775 	*y = x;
776 	return k;
777 	}
778 
779  static Bigint *
780 i2b
781 #ifdef KR_headers
782 	(i) int i;
783 #else
784 	(int i)
785 #endif
786 {
787 	Bigint *b;
788 
789 	b = Balloc(1);
790 	b->x[0] = i;
791 	b->wds = 1;
792 	return b;
793 	}
794 
795  static Bigint *
796 mult
797 #ifdef KR_headers
798 	(a, b) Bigint *a, *b;
799 #else
800 	(Bigint *a, Bigint *b)
801 #endif
802 {
803 	Bigint *c;
804 	int k, wa, wb, wc;
805 	ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
806 	ULong y;
807 #ifdef ULLong
808 	ULLong carry, z;
809 #else
810 	ULong carry, z;
811 #ifdef Pack_32
812 	ULong z2;
813 #endif
814 #endif
815 
816 	if (a->wds < b->wds) {
817 		c = a;
818 		a = b;
819 		b = c;
820 		}
821 	k = a->k;
822 	wa = a->wds;
823 	wb = b->wds;
824 	wc = wa + wb;
825 	if (wc > a->maxwds)
826 		k++;
827 	c = Balloc(k);
828 	for(x = c->x, xa = x + wc; x < xa; x++)
829 		*x = 0;
830 	xa = a->x;
831 	xae = xa + wa;
832 	xb = b->x;
833 	xbe = xb + wb;
834 	xc0 = c->x;
835 #ifdef ULLong
836 	for(; xb < xbe; xc0++) {
837 		if ((y = *xb++)) {
838 			x = xa;
839 			xc = xc0;
840 			carry = 0;
841 			do {
842 				z = *x++ * (ULLong)y + *xc + carry;
843 				carry = z >> 32;
844 				*xc++ = z & FFFFFFFF;
845 				}
846 				while(x < xae);
847 			*xc = carry;
848 			}
849 		}
850 #else
851 #ifdef Pack_32
852 	for(; xb < xbe; xb++, xc0++) {
853 		if (y = *xb & 0xffff) {
854 			x = xa;
855 			xc = xc0;
856 			carry = 0;
857 			do {
858 				z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
859 				carry = z >> 16;
860 				z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
861 				carry = z2 >> 16;
862 				Storeinc(xc, z2, z);
863 				}
864 				while(x < xae);
865 			*xc = carry;
866 			}
867 		if (y = *xb >> 16) {
868 			x = xa;
869 			xc = xc0;
870 			carry = 0;
871 			z2 = *xc;
872 			do {
873 				z = (*x & 0xffff) * y + (*xc >> 16) + carry;
874 				carry = z >> 16;
875 				Storeinc(xc, z, z2);
876 				z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
877 				carry = z2 >> 16;
878 				}
879 				while(x < xae);
880 			*xc = z2;
881 			}
882 		}
883 #else
884 	for(; xb < xbe; xc0++) {
885 		if (y = *xb++) {
886 			x = xa;
887 			xc = xc0;
888 			carry = 0;
889 			do {
890 				z = *x++ * y + *xc + carry;
891 				carry = z >> 16;
892 				*xc++ = z & 0xffff;
893 				}
894 				while(x < xae);
895 			*xc = carry;
896 			}
897 		}
898 #endif
899 #endif
900 	for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
901 	c->wds = wc;
902 	return c;
903 	}
904 
905  static Bigint *p5s;
906 
907  static Bigint *
908 pow5mult
909 #ifdef KR_headers
910 	(b, k) Bigint *b; int k;
911 #else
912 	(Bigint *b, int k)
913 #endif
914 {
915 	Bigint *b1, *p5, *p51;
916 	int i;
917 	static int p05[3] = { 5, 25, 125 };
918 
919 	if ((i = k & 3))
920 		b = multadd(b, p05[i-1], 0);
921 
922 	if (!(k >>= 2))
923 		return b;
924 	if (!(p5 = p5s)) {
925 		/* first time */
926 #ifdef MULTIPLE_THREADS
927 		ACQUIRE_DTOA_LOCK(1);
928 		if (!(p5 = p5s)) {
929 			p5 = p5s = i2b(625);
930 			p5->next = 0;
931 			}
932 		FREE_DTOA_LOCK(1);
933 #else
934 		p5 = p5s = i2b(625);
935 		p5->next = 0;
936 #endif
937 		}
938 	for(;;) {
939 		if (k & 1) {
940 			b1 = mult(b, p5);
941 			Bfree(b);
942 			b = b1;
943 			}
944 		if (!(k >>= 1))
945 			break;
946 		if (!(p51 = p5->next)) {
947 #ifdef MULTIPLE_THREADS
948 			ACQUIRE_DTOA_LOCK(1);
949 			if (!(p51 = p5->next)) {
950 				p51 = p5->next = mult(p5,p5);
951 				p51->next = 0;
952 				}
953 			FREE_DTOA_LOCK(1);
954 #else
955 			p51 = p5->next = mult(p5,p5);
956 			p51->next = 0;
957 #endif
958 			}
959 		p5 = p51;
960 		}
961 	return b;
962 	}
963 
964  static Bigint *
965 lshift
966 #ifdef KR_headers
967 	(b, k) Bigint *b; int k;
968 #else
969 	(Bigint *b, int k)
970 #endif
971 {
972 	int i, k1, n, n1;
973 	Bigint *b1;
974 	ULong *x, *x1, *xe, z;
975 
976 #ifdef Pack_32
977 	n = k >> 5;
978 #else
979 	n = k >> 4;
980 #endif
981 	k1 = b->k;
982 	n1 = n + b->wds + 1;
983 	for(i = b->maxwds; n1 > i; i <<= 1)
984 		k1++;
985 	b1 = Balloc(k1);
986 	x1 = b1->x;
987 	for(i = 0; i < n; i++)
988 		*x1++ = 0;
989 	x = b->x;
990 	xe = x + b->wds;
991 #ifdef Pack_32
992 	if (k &= 0x1f) {
993 		k1 = 32 - k;
994 		z = 0;
995 		do {
996 			*x1++ = *x << k | z;
997 			z = *x++ >> k1;
998 			}
999 			while(x < xe);
1000 		if ((*x1 = z))
1001 			++n1;
1002 		}
1003 #else
1004 	if (k &= 0xf) {
1005 		k1 = 16 - k;
1006 		z = 0;
1007 		do {
1008 			*x1++ = *x << k  & 0xffff | z;
1009 			z = *x++ >> k1;
1010 			}
1011 			while(x < xe);
1012 		if (*x1 = z)
1013 			++n1;
1014 		}
1015 #endif
1016 	else do
1017 		*x1++ = *x++;
1018 		while(x < xe);
1019 	b1->wds = n1 - 1;
1020 	Bfree(b);
1021 	return b1;
1022 	}
1023 
1024  static int
1025 cmp
1026 #ifdef KR_headers
1027 	(a, b) Bigint *a, *b;
1028 #else
1029 	(Bigint *a, Bigint *b)
1030 #endif
1031 {
1032 	ULong *xa, *xa0, *xb, *xb0;
1033 	int i, j;
1034 
1035 	i = a->wds;
1036 	j = b->wds;
1037 #ifdef DEBUG
1038 	if (i > 1 && !a->x[i-1])
1039 		Bug("cmp called with a->x[a->wds-1] == 0");
1040 	if (j > 1 && !b->x[j-1])
1041 		Bug("cmp called with b->x[b->wds-1] == 0");
1042 #endif
1043 	if (i -= j)
1044 		return i;
1045 	xa0 = a->x;
1046 	xa = xa0 + j;
1047 	xb0 = b->x;
1048 	xb = xb0 + j;
1049 	for(;;) {
1050 		if (*--xa != *--xb)
1051 			return *xa < *xb ? -1 : 1;
1052 		if (xa <= xa0)
1053 			break;
1054 		}
1055 	return 0;
1056 	}
1057 
1058  static Bigint *
1059 diff
1060 #ifdef KR_headers
1061 	(a, b) Bigint *a, *b;
1062 #else
1063 	(Bigint *a, Bigint *b)
1064 #endif
1065 {
1066 	Bigint *c;
1067 	int i, wa, wb;
1068 	ULong *xa, *xae, *xb, *xbe, *xc;
1069 #ifdef ULLong
1070 	ULLong borrow, y;
1071 #else
1072 	ULong borrow, y;
1073 #ifdef Pack_32
1074 	ULong z;
1075 #endif
1076 #endif
1077 
1078 	i = cmp(a,b);
1079 	if (!i) {
1080 		c = Balloc(0);
1081 		c->wds = 1;
1082 		c->x[0] = 0;
1083 		return c;
1084 		}
1085 	if (i < 0) {
1086 		c = a;
1087 		a = b;
1088 		b = c;
1089 		i = 1;
1090 		}
1091 	else
1092 		i = 0;
1093 	c = Balloc(a->k);
1094 	c->sign = i;
1095 	wa = a->wds;
1096 	xa = a->x;
1097 	xae = xa + wa;
1098 	wb = b->wds;
1099 	xb = b->x;
1100 	xbe = xb + wb;
1101 	xc = c->x;
1102 	borrow = 0;
1103 #ifdef ULLong
1104 	do {
1105 		y = (ULLong)*xa++ - *xb++ - borrow;
1106 		borrow = y >> 32 & (ULong)1;
1107 		*xc++ = y & FFFFFFFF;
1108 		}
1109 		while(xb < xbe);
1110 	while(xa < xae) {
1111 		y = *xa++ - borrow;
1112 		borrow = y >> 32 & (ULong)1;
1113 		*xc++ = y & FFFFFFFF;
1114 		}
1115 #else
1116 #ifdef Pack_32
1117 	do {
1118 		y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
1119 		borrow = (y & 0x10000) >> 16;
1120 		z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
1121 		borrow = (z & 0x10000) >> 16;
1122 		Storeinc(xc, z, y);
1123 		}
1124 		while(xb < xbe);
1125 	while(xa < xae) {
1126 		y = (*xa & 0xffff) - borrow;
1127 		borrow = (y & 0x10000) >> 16;
1128 		z = (*xa++ >> 16) - borrow;
1129 		borrow = (z & 0x10000) >> 16;
1130 		Storeinc(xc, z, y);
1131 		}
1132 #else
1133 	do {
1134 		y = *xa++ - *xb++ - borrow;
1135 		borrow = (y & 0x10000) >> 16;
1136 		*xc++ = y & 0xffff;
1137 		}
1138 		while(xb < xbe);
1139 	while(xa < xae) {
1140 		y = *xa++ - borrow;
1141 		borrow = (y & 0x10000) >> 16;
1142 		*xc++ = y & 0xffff;
1143 		}
1144 #endif
1145 #endif
1146 	while(!*--xc)
1147 		wa--;
1148 	c->wds = wa;
1149 	return c;
1150 	}
1151 
1152  static double
1153 ulp
1154 #ifdef KR_headers
1155 	(x) U *x;
1156 #else
1157 	(U *x)
1158 #endif
1159 {
1160 	Long L;
1161 	U u;
1162 
1163 	L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
1164 #ifndef Avoid_Underflow
1165 #ifndef Sudden_Underflow
1166 	if (L > 0) {
1167 #endif
1168 #endif
1169 #ifdef IBM
1170 		L |= Exp_msk1 >> 4;
1171 #endif
1172 		word0(&u) = L;
1173 		word1(&u) = 0;
1174 #ifndef Avoid_Underflow
1175 #ifndef Sudden_Underflow
1176 		}
1177 	else {
1178 		L = -L >> Exp_shift;
1179 		if (L < Exp_shift) {
1180 			word0(&u) = 0x80000 >> L;
1181 			word1(&u) = 0;
1182 			}
1183 		else {
1184 			word0(&u) = 0;
1185 			L -= Exp_shift;
1186 			word1(&u) = L >= 31 ? 1 : 1 << 31 - L;
1187 			}
1188 		}
1189 #endif
1190 #endif
1191 	return dval(&u);
1192 	}
1193 
1194  static double
1195 b2d
1196 #ifdef KR_headers
1197 	(a, e) Bigint *a; int *e;
1198 #else
1199 	(Bigint *a, int *e)
1200 #endif
1201 {
1202 	ULong *xa, *xa0, w, y, z;
1203 	int k;
1204 	U d;
1205 #ifdef VAX
1206 	ULong d0, d1;
1207 #else
1208 #define d0 word0(&d)
1209 #define d1 word1(&d)
1210 #endif
1211 
1212 	xa0 = a->x;
1213 	xa = xa0 + a->wds;
1214 	y = *--xa;
1215 #ifdef DEBUG
1216 	if (!y) Bug("zero y in b2d");
1217 #endif
1218 	k = hi0bits(y);
1219 	*e = 32 - k;
1220 #ifdef Pack_32
1221 	if (k < Ebits) {
1222 		d0 = Exp_1 | y >> (Ebits - k);
1223 		w = xa > xa0 ? *--xa : 0;
1224 		d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
1225 		goto ret_d;
1226 		}
1227 	z = xa > xa0 ? *--xa : 0;
1228 	if (k -= Ebits) {
1229 		d0 = Exp_1 | y << k | z >> (32 - k);
1230 		y = xa > xa0 ? *--xa : 0;
1231 		d1 = z << k | y >> (32 - k);
1232 		}
1233 	else {
1234 		d0 = Exp_1 | y;
1235 		d1 = z;
1236 		}
1237 #else
1238 	if (k < Ebits + 16) {
1239 		z = xa > xa0 ? *--xa : 0;
1240 		d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
1241 		w = xa > xa0 ? *--xa : 0;
1242 		y = xa > xa0 ? *--xa : 0;
1243 		d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
1244 		goto ret_d;
1245 		}
1246 	z = xa > xa0 ? *--xa : 0;
1247 	w = xa > xa0 ? *--xa : 0;
1248 	k -= Ebits + 16;
1249 	d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
1250 	y = xa > xa0 ? *--xa : 0;
1251 	d1 = w << k + 16 | y << k;
1252 #endif
1253  ret_d:
1254 #ifdef VAX
1255 	word0(&d) = d0 >> 16 | d0 << 16;
1256 	word1(&d) = d1 >> 16 | d1 << 16;
1257 #else
1258 #undef d0
1259 #undef d1
1260 #endif
1261 	return dval(&d);
1262 	}
1263 
1264  static Bigint *
1265 d2b
1266 #ifdef KR_headers
1267 	(d, e, bits) U *d; int *e, *bits;
1268 #else
1269 	(U *d, int *e, int *bits)
1270 #endif
1271 {
1272 	Bigint *b;
1273 	int de, k;
1274 	ULong *x, y, z;
1275 #ifndef Sudden_Underflow
1276 	int i;
1277 #endif
1278 #ifdef VAX
1279 	ULong d0, d1;
1280 	d0 = word0(d) >> 16 | word0(d) << 16;
1281 	d1 = word1(d) >> 16 | word1(d) << 16;
1282 #else
1283 #define d0 word0(d)
1284 #define d1 word1(d)
1285 #endif
1286 
1287 #ifdef Pack_32
1288 	b = Balloc(1);
1289 #else
1290 	b = Balloc(2);
1291 #endif
1292 	x = b->x;
1293 
1294 	z = d0 & Frac_mask;
1295 	d0 &= 0x7fffffff;	/* clear sign bit, which we ignore */
1296 #ifdef Sudden_Underflow
1297 	de = (int)(d0 >> Exp_shift);
1298 #ifndef IBM
1299 	z |= Exp_msk11;
1300 #endif
1301 #else
1302 	if ((de = (int)(d0 >> Exp_shift)))
1303 		z |= Exp_msk1;
1304 #endif
1305 #ifdef Pack_32
1306 	if ((y = d1)) {
1307 		if ((k = lo0bits(&y))) {
1308 			x[0] = y | z << (32 - k);
1309 			z >>= k;
1310 			}
1311 		else
1312 			x[0] = y;
1313 #ifndef Sudden_Underflow
1314 		i =
1315 #endif
1316 		    b->wds = (x[1] = z) ? 2 : 1;
1317 		}
1318 	else {
1319 		k = lo0bits(&z);
1320 		x[0] = z;
1321 #ifndef Sudden_Underflow
1322 		i =
1323 #endif
1324 		    b->wds = 1;
1325 		k += 32;
1326 		}
1327 #else
1328 	if (y = d1) {
1329 		if (k = lo0bits(&y))
1330 			if (k >= 16) {
1331 				x[0] = y | z << 32 - k & 0xffff;
1332 				x[1] = z >> k - 16 & 0xffff;
1333 				x[2] = z >> k;
1334 				i = 2;
1335 				}
1336 			else {
1337 				x[0] = y & 0xffff;
1338 				x[1] = y >> 16 | z << 16 - k & 0xffff;
1339 				x[2] = z >> k & 0xffff;
1340 				x[3] = z >> k+16;
1341 				i = 3;
1342 				}
1343 		else {
1344 			x[0] = y & 0xffff;
1345 			x[1] = y >> 16;
1346 			x[2] = z & 0xffff;
1347 			x[3] = z >> 16;
1348 			i = 3;
1349 			}
1350 		}
1351 	else {
1352 #ifdef DEBUG
1353 		if (!z)
1354 			Bug("Zero passed to d2b");
1355 #endif
1356 		k = lo0bits(&z);
1357 		if (k >= 16) {
1358 			x[0] = z;
1359 			i = 0;
1360 			}
1361 		else {
1362 			x[0] = z & 0xffff;
1363 			x[1] = z >> 16;
1364 			i = 1;
1365 			}
1366 		k += 32;
1367 		}
1368 	while(!x[i])
1369 		--i;
1370 	b->wds = i + 1;
1371 #endif
1372 #ifndef Sudden_Underflow
1373 	if (de) {
1374 #endif
1375 #ifdef IBM
1376 		*e = (de - Bias - (P-1) << 2) + k;
1377 		*bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
1378 #else
1379 		*e = de - Bias - (P-1) + k;
1380 		*bits = P - k;
1381 #endif
1382 #ifndef Sudden_Underflow
1383 		}
1384 	else {
1385 		*e = de - Bias - (P-1) + 1 + k;
1386 #ifdef Pack_32
1387 		*bits = 32*i - hi0bits(x[i-1]);
1388 #else
1389 		*bits = (i+2)*16 - hi0bits(x[i]);
1390 #endif
1391 		}
1392 #endif
1393 	return b;
1394 	}
1395 #undef d0
1396 #undef d1
1397 
1398  static double
1399 ratio
1400 #ifdef KR_headers
1401 	(a, b) Bigint *a, *b;
1402 #else
1403 	(Bigint *a, Bigint *b)
1404 #endif
1405 {
1406 	U da, db;
1407 	int k, ka, kb;
1408 
1409 	dval(&da) = b2d(a, &ka);
1410 	dval(&db) = b2d(b, &kb);
1411 #ifdef Pack_32
1412 	k = ka - kb + 32*(a->wds - b->wds);
1413 #else
1414 	k = ka - kb + 16*(a->wds - b->wds);
1415 #endif
1416 #ifdef IBM
1417 	if (k > 0) {
1418 		word0(&da) += (k >> 2)*Exp_msk1;
1419 		if (k &= 3)
1420 			dval(&da) *= 1 << k;
1421 		}
1422 	else {
1423 		k = -k;
1424 		word0(&db) += (k >> 2)*Exp_msk1;
1425 		if (k &= 3)
1426 			dval(&db) *= 1 << k;
1427 		}
1428 #else
1429 	if (k > 0)
1430 		word0(&da) += k*Exp_msk1;
1431 	else {
1432 		k = -k;
1433 		word0(&db) += k*Exp_msk1;
1434 		}
1435 #endif
1436 	return dval(&da) / dval(&db);
1437 	}
1438 
1439  static CONST double
1440 tens[] = {
1441 		1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1442 		1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1443 		1e20, 1e21, 1e22
1444 #ifdef VAX
1445 		, 1e23, 1e24
1446 #endif
1447 		};
1448 
1449  static CONST double
1450 #ifdef IEEE_Arith
1451 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1452 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
1453 #ifdef Avoid_Underflow
1454 		9007199254740992.*9007199254740992.e-256
1455 		/* = 2^106 * 1e-256 */
1456 #else
1457 		1e-256
1458 #endif
1459 		};
1460 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
1461 /* flag unnecessarily.  It leads to a song and dance at the end of strtod. */
1462 #define Scale_Bit 0x10
1463 #define n_bigtens 5
1464 #else
1465 #ifdef IBM
1466 bigtens[] = { 1e16, 1e32, 1e64 };
1467 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
1468 #define n_bigtens 3
1469 #else
1470 bigtens[] = { 1e16, 1e32 };
1471 static CONST double tinytens[] = { 1e-16, 1e-32 };
1472 #define n_bigtens 2
1473 #endif
1474 #endif
1475 
1476 #undef Need_Hexdig
1477 #ifdef INFNAN_CHECK
1478 #ifndef No_Hex_NaN
1479 #define Need_Hexdig
1480 #endif
1481 #endif
1482 
1483 #ifndef Need_Hexdig
1484 #ifndef NO_HEX_FP
1485 #define Need_Hexdig
1486 #endif
1487 #endif
1488 
1489 #ifdef Need_Hexdig /*{*/
1490 static unsigned char hexdig[256];
1491 
1492  static void
1493 #ifdef KR_headers
1494 htinit(h, s, inc) unsigned char *h; unsigned char *s; int inc;
1495 #else
1496 htinit(unsigned char *h, unsigned char *s, int inc)
1497 #endif
1498 {
1499 	int i, j;
1500 	for(i = 0; (j = s[i]) !=0; i++)
1501 		h[j] = i + inc;
1502 	}
1503 
1504  static void
1505 #ifdef KR_headers
hexdig_init()1506 hexdig_init()
1507 #else
1508 hexdig_init(void)
1509 #endif
1510 {
1511 #define USC (unsigned char *)
1512 	htinit(hexdig, USC "0123456789", 0x10);
1513 	htinit(hexdig, USC "abcdef", 0x10 + 10);
1514 	htinit(hexdig, USC "ABCDEF", 0x10 + 10);
1515 	}
1516 #endif /* } Need_Hexdig */
1517 
1518 #ifdef INFNAN_CHECK
1519 
1520 #ifndef NAN_WORD0
1521 #define NAN_WORD0 0x7ff80000
1522 #endif
1523 
1524 #ifndef NAN_WORD1
1525 #define NAN_WORD1 0
1526 #endif
1527 
1528  static int
1529 match
1530 #ifdef KR_headers
1531 	(sp, t) char **sp, *t;
1532 #else
1533 	(CONST char **sp, CONST char *t)
1534 #endif
1535 {
1536 	int c, d;
1537 	CONST char *s = *sp;
1538 
1539 	while((d = *t++)) {
1540 		if ((c = *++s) >= 'A' && c <= 'Z')
1541 			c += 'a' - 'A';
1542 		if (c != d)
1543 			return 0;
1544 		}
1545 	*sp = s + 1;
1546 	return 1;
1547 	}
1548 
1549 #ifndef No_Hex_NaN
1550  static void
1551 hexnan
1552 #ifdef KR_headers
1553 	(rvp, sp) U *rvp; CONST char **sp;
1554 #else
1555 	(U *rvp, CONST char **sp)
1556 #endif
1557 {
1558 	ULong c, x[2];
1559 	CONST char *s;
1560 	int c1, havedig, udx0, xshift;
1561 
1562 	if (!hexdig['0'])
1563 		hexdig_init();
1564 	x[0] = x[1] = 0;
1565 	havedig = xshift = 0;
1566 	udx0 = 1;
1567 	s = *sp;
1568 	/* allow optional initial 0x or 0X */
1569 	while((c = *(CONST unsigned char*)(s+1)) && c <= ' ')
1570 		++s;
1571 	if (s[1] == '0' && (s[2] == 'x' || s[2] == 'X'))
1572 		s += 2;
1573 	while((c = *(CONST unsigned char*)++s)) {
1574 		if ((c1 = hexdig[c]))
1575 			c  = c1 & 0xf;
1576 		else if (c <= ' ') {
1577 			if (udx0 && havedig) {
1578 				udx0 = 0;
1579 				xshift = 1;
1580 				}
1581 			continue;
1582 			}
1583 #ifdef GDTOA_NON_PEDANTIC_NANCHECK
1584 		else if (/*(*/ c == ')' && havedig) {
1585 			*sp = s + 1;
1586 			break;
1587 			}
1588 		else
1589 			return;	/* invalid form: don't change *sp */
1590 #else
1591 		else {
1592 			do {
1593 				if (/*(*/ c == ')') {
1594 					*sp = s + 1;
1595 					break;
1596 					}
1597 				} while((c = *++s));
1598 			break;
1599 			}
1600 #endif
1601 		havedig = 1;
1602 		if (xshift) {
1603 			xshift = 0;
1604 			x[0] = x[1];
1605 			x[1] = 0;
1606 			}
1607 		if (udx0)
1608 			x[0] = (x[0] << 4) | (x[1] >> 28);
1609 		x[1] = (x[1] << 4) | c;
1610 		}
1611 	if ((x[0] &= 0xfffff) || x[1]) {
1612 		word0(rvp) = Exp_mask | x[0];
1613 		word1(rvp) = x[1];
1614 		}
1615 	}
1616 #endif /*No_Hex_NaN*/
1617 #endif /* INFNAN_CHECK */
1618 
1619 #ifdef Pack_32
1620 #define ULbits 32
1621 #define kshift 5
1622 #define kmask 31
1623 #else
1624 #define ULbits 16
1625 #define kshift 4
1626 #define kmask 15
1627 #endif
1628 #ifndef NO_HEX_FP /*{*/
1629 
1630  static void
1631 #ifdef KR_headers
1632 rshift(b, k) Bigint *b; int k;
1633 #else
1634 rshift(Bigint *b, int k)
1635 #endif
1636 {
1637 	ULong *x, *x1, *xe, y;
1638 	int n;
1639 
1640 	x = x1 = b->x;
1641 	n = k >> kshift;
1642 	if (n < b->wds) {
1643 		xe = x + b->wds;
1644 		x += n;
1645 		if (k &= kmask) {
1646 			n = 32 - k;
1647 			y = *x++ >> k;
1648 			while(x < xe) {
1649 				*x1++ = (y | (*x << n)) & 0xffffffff;
1650 				y = *x++ >> k;
1651 				}
1652 			if ((*x1 = y) !=0)
1653 				x1++;
1654 			}
1655 		else
1656 			while(x < xe)
1657 				*x1++ = *x++;
1658 		}
1659 	if ((b->wds = x1 - b->x) == 0)
1660 		b->x[0] = 0;
1661 	}
1662 
1663  static ULong
1664 #ifdef KR_headers
1665 any_on(b, k) Bigint *b; int k;
1666 #else
1667 any_on(Bigint *b, int k)
1668 #endif
1669 {
1670 	int n, nwds;
1671 	ULong *x, *x0, x1, x2;
1672 
1673 	x = b->x;
1674 	nwds = b->wds;
1675 	n = k >> kshift;
1676 	if (n > nwds)
1677 		n = nwds;
1678 	else if (n < nwds && (k &= kmask)) {
1679 		x1 = x2 = x[n];
1680 		x1 >>= k;
1681 		x1 <<= k;
1682 		if (x1 != x2)
1683 			return 1;
1684 		}
1685 	x0 = x;
1686 	x += n;
1687 	while(x > x0)
1688 		if (*--x)
1689 			return 1;
1690 	return 0;
1691 	}
1692 
1693 enum {	/* rounding values: same as FLT_ROUNDS */
1694 	Round_zero = 0,
1695 	Round_near = 1,
1696 	Round_up = 2,
1697 	Round_down = 3
1698 	};
1699 
1700  static Bigint *
1701 #ifdef KR_headers
1702 increment(b) Bigint *b;
1703 #else
1704 increment(Bigint *b)
1705 #endif
1706 {
1707 	ULong *x, *xe;
1708 	Bigint *b1;
1709 
1710 	x = b->x;
1711 	xe = x + b->wds;
1712 	do {
1713 		if (*x < (ULong)0xffffffffL) {
1714 			++*x;
1715 			return b;
1716 			}
1717 		*x++ = 0;
1718 		} while(x < xe);
1719 	{
1720 		if (b->wds >= b->maxwds) {
1721 			b1 = Balloc(b->k+1);
1722 			Bcopy(b1,b);
1723 			Bfree(b);
1724 			b = b1;
1725 			}
1726 		b->x[b->wds++] = 1;
1727 		}
1728 	return b;
1729 	}
1730 
1731  void
1732 #ifdef KR_headers
1733 gethex(sp, rvp, rounding, sign)
1734 	CONST char **sp; U *rvp; int rounding, sign;
1735 #else
1736 gethex( CONST char **sp, U *rvp, int rounding, int sign)
1737 #endif
1738 {
1739 	Bigint *b;
1740 	CONST unsigned char *decpt, *s0, *s, *s1;
1741 	Long e, e1;
1742 	ULong L, lostbits, *x;
1743 	int big, denorm, esign, havedig, k, n, nbits, up, zret;
1744 #ifdef IBM
1745 	int j;
1746 #endif
1747 	enum {
1748 #ifdef IEEE_Arith /*{{*/
1749 		emax = 0x7fe - Bias - P + 1,
1750 		emin = Emin - P + 1
1751 #else /*}{*/
1752 		emin = Emin - P,
1753 #ifdef VAX
1754 		emax = 0x7ff - Bias - P + 1
1755 #endif
1756 #ifdef IBM
1757 		emax = 0x7f - Bias - P
1758 #endif
1759 #endif /*}}*/
1760 		};
1761 #ifdef USE_LOCALE
1762 	int i;
1763 #ifdef NO_LOCALE_CACHE
1764 	const unsigned char *decimalpoint = (unsigned char*)
1765 		localeconv()->decimal_point;
1766 #else
1767 	const unsigned char *decimalpoint;
1768 	static unsigned char *decimalpoint_cache;
1769 	if (!(s0 = decimalpoint_cache)) {
1770 		s0 = (unsigned char*)localeconv()->decimal_point;
1771 		if ((decimalpoint_cache = (unsigned char*)
1772 				MALLOC(strlen((CONST char*)s0) + 1))) {
1773 			strcpy((char*)decimalpoint_cache, (CONST char*)s0);
1774 			s0 = decimalpoint_cache;
1775 			}
1776 		}
1777 	decimalpoint = s0;
1778 #endif
1779 #endif
1780 
1781 	if (!hexdig['0'])
1782 		hexdig_init();
1783 	havedig = 0;
1784 	s0 = *(CONST unsigned char **)sp + 2;
1785 	while(s0[havedig] == '0')
1786 		havedig++;
1787 	s0 += havedig;
1788 	s = s0;
1789 	decpt = 0;
1790 	zret = 0;
1791 	e = 0;
1792 	if (hexdig[*s])
1793 		havedig++;
1794 	else {
1795 		zret = 1;
1796 #ifdef USE_LOCALE
1797 		for(i = 0; decimalpoint[i]; ++i) {
1798 			if (s[i] != decimalpoint[i])
1799 				goto pcheck;
1800 			}
1801 		decpt = s += i;
1802 #else
1803 		if (*s != '.')
1804 			goto pcheck;
1805 		decpt = ++s;
1806 #endif
1807 		if (!hexdig[*s])
1808 			goto pcheck;
1809 		while(*s == '0')
1810 			s++;
1811 		if (hexdig[*s])
1812 			zret = 0;
1813 		havedig = 1;
1814 		s0 = s;
1815 		}
1816 	while(hexdig[*s])
1817 		s++;
1818 #ifdef USE_LOCALE
1819 	if (*s == *decimalpoint && !decpt) {
1820 		for(i = 1; decimalpoint[i]; ++i) {
1821 			if (s[i] != decimalpoint[i])
1822 				goto pcheck;
1823 			}
1824 		decpt = s += i;
1825 #else
1826 	if (*s == '.' && !decpt) {
1827 		decpt = ++s;
1828 #endif
1829 		while(hexdig[*s])
1830 			s++;
1831 		}/*}*/
1832 	if (decpt)
1833 		e = -(((Long)(s-decpt)) << 2);
1834  pcheck:
1835 	s1 = s;
1836 	big = esign = 0;
1837 	switch(*s) {
1838 	  case 'p':
1839 	  case 'P':
1840 		switch(*++s) {
1841 		  case '-':
1842 			esign = 1;
1843 			/* no break */
1844 		  case '+':
1845 			s++;
1846 		  }
1847 		if ((n = hexdig[*s]) == 0 || n > 0x19) {
1848 			s = s1;
1849 			break;
1850 			}
1851 		e1 = n - 0x10;
1852 		while((n = hexdig[*++s]) !=0 && n <= 0x19) {
1853 			if (e1 & 0xf8000000)
1854 				big = 1;
1855 			e1 = 10*e1 + n - 0x10;
1856 			}
1857 		if (esign)
1858 			e1 = -e1;
1859 		e += e1;
1860 	  }
1861 	*sp = (char*)s;
1862 	if (!havedig)
1863 		*sp = (char*)s0 - 1;
1864 	if (zret)
1865 		goto retz1;
1866 	if (big) {
1867 		if (esign) {
1868 #ifdef IEEE_Arith
1869 			switch(rounding) {
1870 			  case Round_up:
1871 				if (sign)
1872 					break;
1873 				goto ret_tiny;
1874 			  case Round_down:
1875 				if (!sign)
1876 					break;
1877 				goto ret_tiny;
1878 			  }
1879 #endif
1880 			goto retz;
1881 #ifdef IEEE_Arith
1882  ret_tiny:
1883 #ifndef NO_ERRNO
1884 			errno = ERANGE;
1885 #endif
1886 			word0(rvp) = 0;
1887 			word1(rvp) = 1;
1888 			return;
1889 #endif /* IEEE_Arith */
1890 			}
1891 		switch(rounding) {
1892 		  case Round_near:
1893 			goto ovfl1;
1894 		  case Round_up:
1895 			if (!sign)
1896 				goto ovfl1;
1897 			goto ret_big;
1898 		  case Round_down:
1899 			if (sign)
1900 				goto ovfl1;
1901 			goto ret_big;
1902 		  }
1903  ret_big:
1904 		word0(rvp) = Big0;
1905 		word1(rvp) = Big1;
1906 		return;
1907 		}
1908 	n = s1 - s0 - 1;
1909 	for(k = 0; n > (1 << (kshift-2)) - 1; n >>= 1)
1910 		k++;
1911 	b = Balloc(k);
1912 	x = b->x;
1913 	n = 0;
1914 	L = 0;
1915 #ifdef USE_LOCALE
1916 	for(i = 0; decimalpoint[i+1]; ++i);
1917 #endif
1918 	while(s1 > s0) {
1919 #ifdef USE_LOCALE
1920 		if (*--s1 == decimalpoint[i]) {
1921 			s1 -= i;
1922 			continue;
1923 			}
1924 #else
1925 		if (*--s1 == '.')
1926 			continue;
1927 #endif
1928 		if (n == ULbits) {
1929 			*x++ = L;
1930 			L = 0;
1931 			n = 0;
1932 			}
1933 		L |= (hexdig[*s1] & 0x0f) << n;
1934 		n += 4;
1935 		}
1936 	*x++ = L;
1937 	b->wds = n = x - b->x;
1938 	n = ULbits*n - hi0bits(L);
1939 	nbits = Nbits;
1940 	lostbits = 0;
1941 	x = b->x;
1942 	if (n > nbits) {
1943 		n -= nbits;
1944 		if (any_on(b,n)) {
1945 			lostbits = 1;
1946 			k = n - 1;
1947 			if (x[k>>kshift] & 1 << (k & kmask)) {
1948 				lostbits = 2;
1949 				if (k > 0 && any_on(b,k))
1950 					lostbits = 3;
1951 				}
1952 			}
1953 		rshift(b, n);
1954 		e += n;
1955 		}
1956 	else if (n < nbits) {
1957 		n = nbits - n;
1958 		b = lshift(b, n);
1959 		e -= n;
1960 		x = b->x;
1961 		}
1962 	if (e > Emax) {
1963  ovfl:
1964 		Bfree(b);
1965  ovfl1:
1966 #ifndef NO_ERRNO
1967 		errno = ERANGE;
1968 #endif
1969 		word0(rvp) = Exp_mask;
1970 		word1(rvp) = 0;
1971 		return;
1972 		}
1973 	denorm = 0;
1974 	if (e < emin) {
1975 		denorm = 1;
1976 		n = emin - e;
1977 		if (n >= nbits) {
1978 #ifdef IEEE_Arith /*{*/
1979 			switch (rounding) {
1980 			  case Round_near:
1981 				if (n == nbits && (n < 2 || any_on(b,n-1)))
1982 					goto ret_tiny;
1983 				break;
1984 			  case Round_up:
1985 				if (!sign)
1986 					goto ret_tiny;
1987 				break;
1988 			  case Round_down:
1989 				if (sign)
1990 					goto ret_tiny;
1991 			  }
1992 #endif /* } IEEE_Arith */
1993 			Bfree(b);
1994  retz:
1995 #ifndef NO_ERRNO
1996 			errno = ERANGE;
1997 #endif
1998  retz1:
1999 			rvp->d = 0.;
2000 			return;
2001 			}
2002 		k = n - 1;
2003 		if (lostbits)
2004 			lostbits = 1;
2005 		else if (k > 0)
2006 			lostbits = any_on(b,k);
2007 		if (x[k>>kshift] & 1 << (k & kmask))
2008 			lostbits |= 2;
2009 		nbits -= n;
2010 		rshift(b,n);
2011 		e = emin;
2012 		}
2013 	if (lostbits) {
2014 		up = 0;
2015 		switch(rounding) {
2016 		  case Round_zero:
2017 			break;
2018 		  case Round_near:
2019 			if (lostbits & 2
2020 			 && (lostbits & 1) | (x[0] & 1))
2021 				up = 1;
2022 			break;
2023 		  case Round_up:
2024 			up = 1 - sign;
2025 			break;
2026 		  case Round_down:
2027 			up = sign;
2028 		  }
2029 		if (up) {
2030 			k = b->wds;
2031 			b = increment(b);
2032 			x = b->x;
2033 			if (denorm) {
2034 #if 0
2035 				if (nbits == Nbits - 1
2036 				 && x[nbits >> kshift] & 1 << (nbits & kmask))
2037 					denorm = 0; /* not currently used */
2038 #endif
2039 				}
2040 			else if (b->wds > k
2041 			 || ((n = nbits & kmask) !=0
2042 			     && hi0bits(x[k-1]) < 32-n)) {
2043 				rshift(b,1);
2044 				if (++e > Emax)
2045 					goto ovfl;
2046 				}
2047 			}
2048 		}
2049 #ifdef IEEE_Arith
2050 	if (denorm)
2051 		word0(rvp) = b->wds > 1 ? b->x[1] & ~0x100000 : 0;
2052 	else
2053 		word0(rvp) = (b->x[1] & ~0x100000) | ((e + 0x3ff + 52) << 20);
2054 	word1(rvp) = b->x[0];
2055 #endif
2056 #ifdef IBM
2057 	if ((j = e & 3)) {
2058 		k = b->x[0] & ((1 << j) - 1);
2059 		rshift(b,j);
2060 		if (k) {
2061 			switch(rounding) {
2062 			  case Round_up:
2063 				if (!sign)
2064 					increment(b);
2065 				break;
2066 			  case Round_down:
2067 				if (sign)
2068 					increment(b);
2069 				break;
2070 			  case Round_near:
2071 				j = 1 << (j-1);
2072 				if (k & j && ((k & (j-1)) | lostbits))
2073 					increment(b);
2074 			  }
2075 			}
2076 		}
2077 	e >>= 2;
2078 	word0(rvp) = b->x[1] | ((e + 65 + 13) << 24);
2079 	word1(rvp) = b->x[0];
2080 #endif
2081 #ifdef VAX
2082 	/* The next two lines ignore swap of low- and high-order 2 bytes. */
2083 	/* word0(rvp) = (b->x[1] & ~0x800000) | ((e + 129 + 55) << 23); */
2084 	/* word1(rvp) = b->x[0]; */
2085 	word0(rvp) = ((b->x[1] & ~0x800000) >> 16) | ((e + 129 + 55) << 7) | (b->x[1] << 16);
2086 	word1(rvp) = (b->x[0] >> 16) | (b->x[0] << 16);
2087 #endif
2088 	Bfree(b);
2089 	}
2090 #endif /*}!NO_HEX_FP*/
2091 
2092  static int
2093 #ifdef KR_headers
2094 dshift(b, p2) Bigint *b; int p2;
2095 #else
2096 dshift(Bigint *b, int p2)
2097 #endif
2098 {
2099 	int rv = hi0bits(b->x[b->wds-1]) - 4;
2100 	if (p2 > 0)
2101 		rv -= p2;
2102 	return rv & kmask;
2103 	}
2104 
2105  static int
2106 quorem
2107 #ifdef KR_headers
2108 	(b, S) Bigint *b, *S;
2109 #else
2110 	(Bigint *b, Bigint *S)
2111 #endif
2112 {
2113 	int n;
2114 	ULong *bx, *bxe, q, *sx, *sxe;
2115 #ifdef ULLong
2116 	ULLong borrow, carry, y, ys;
2117 #else
2118 	ULong borrow, carry, y, ys;
2119 #ifdef Pack_32
2120 	ULong si, z, zs;
2121 #endif
2122 #endif
2123 
2124 	n = S->wds;
2125 #ifdef DEBUG
2126 	/*debug*/ if (b->wds > n)
2127 	/*debug*/	Bug("oversize b in quorem");
2128 #endif
2129 	if (b->wds < n)
2130 		return 0;
2131 	sx = S->x;
2132 	sxe = sx + --n;
2133 	bx = b->x;
2134 	bxe = bx + n;
2135 	q = *bxe / (*sxe + 1);	/* ensure q <= true quotient */
2136 #ifdef DEBUG
2137 	/*debug*/ if (q > 9)
2138 	/*debug*/	Bug("oversized quotient in quorem");
2139 #endif
2140 	if (q) {
2141 		borrow = 0;
2142 		carry = 0;
2143 		do {
2144 #ifdef ULLong
2145 			ys = *sx++ * (ULLong)q + carry;
2146 			carry = ys >> 32;
2147 			y = *bx - (ys & FFFFFFFF) - borrow;
2148 			borrow = y >> 32 & (ULong)1;
2149 			*bx++ = y & FFFFFFFF;
2150 #else
2151 #ifdef Pack_32
2152 			si = *sx++;
2153 			ys = (si & 0xffff) * q + carry;
2154 			zs = (si >> 16) * q + (ys >> 16);
2155 			carry = zs >> 16;
2156 			y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
2157 			borrow = (y & 0x10000) >> 16;
2158 			z = (*bx >> 16) - (zs & 0xffff) - borrow;
2159 			borrow = (z & 0x10000) >> 16;
2160 			Storeinc(bx, z, y);
2161 #else
2162 			ys = *sx++ * q + carry;
2163 			carry = ys >> 16;
2164 			y = *bx - (ys & 0xffff) - borrow;
2165 			borrow = (y & 0x10000) >> 16;
2166 			*bx++ = y & 0xffff;
2167 #endif
2168 #endif
2169 			}
2170 			while(sx <= sxe);
2171 		if (!*bxe) {
2172 			bx = b->x;
2173 			while(--bxe > bx && !*bxe)
2174 				--n;
2175 			b->wds = n;
2176 			}
2177 		}
2178 	if (cmp(b, S) >= 0) {
2179 		q++;
2180 		borrow = 0;
2181 		carry = 0;
2182 		bx = b->x;
2183 		sx = S->x;
2184 		do {
2185 #ifdef ULLong
2186 			ys = *sx++ + carry;
2187 			carry = ys >> 32;
2188 			y = *bx - (ys & FFFFFFFF) - borrow;
2189 			borrow = y >> 32 & (ULong)1;
2190 			*bx++ = y & FFFFFFFF;
2191 #else
2192 #ifdef Pack_32
2193 			si = *sx++;
2194 			ys = (si & 0xffff) + carry;
2195 			zs = (si >> 16) + (ys >> 16);
2196 			carry = zs >> 16;
2197 			y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
2198 			borrow = (y & 0x10000) >> 16;
2199 			z = (*bx >> 16) - (zs & 0xffff) - borrow;
2200 			borrow = (z & 0x10000) >> 16;
2201 			Storeinc(bx, z, y);
2202 #else
2203 			ys = *sx++ + carry;
2204 			carry = ys >> 16;
2205 			y = *bx - (ys & 0xffff) - borrow;
2206 			borrow = (y & 0x10000) >> 16;
2207 			*bx++ = y & 0xffff;
2208 #endif
2209 #endif
2210 			}
2211 			while(sx <= sxe);
2212 		bx = b->x;
2213 		bxe = bx + n;
2214 		if (!*bxe) {
2215 			while(--bxe > bx && !*bxe)
2216 				--n;
2217 			b->wds = n;
2218 			}
2219 		}
2220 	return q;
2221 	}
2222 
2223 #ifndef NO_STRTOD_BIGCOMP
2224 
2225  static void
2226 bigcomp
2227 #ifdef KR_headers
2228 	(rv, s0, bc)
2229 	U *rv; CONST char *s0; BCinfo *bc;
2230 #else
2231 	(U *rv, CONST char *s0, BCinfo *bc)
2232 #endif
2233 {
2234 	Bigint *b, *d;
2235 	int b2, bbits, d2, dd, dig, dsign, i, j, nd, nd0, p2, p5, speccase;
2236 
2237 	dsign = bc->dsign;
2238 	nd = bc->nd;
2239 	nd0 = bc->nd0;
2240 	p5 = nd + bc->e0 - 1;
2241 	dd = speccase = 0;
2242 #ifndef Sudden_Underflow
2243 	if (rv->d == 0.) {	/* special case: value near underflow-to-zero */
2244 				/* threshold was rounded to zero */
2245 		b = i2b(1);
2246 		p2 = Emin - P + 1;
2247 		bbits = 1;
2248 #ifdef Avoid_Underflow
2249 		word0(rv) = (P+2) << Exp_shift;
2250 #else
2251 		word1(rv) = 1;
2252 #endif
2253 		i = 0;
2254 #ifdef Honor_FLT_ROUNDS
2255 		if (bc->rounding == 1)
2256 #endif
2257 			{
2258 			speccase = 1;
2259 			--p2;
2260 			dsign = 0;
2261 			goto have_i;
2262 			}
2263 		}
2264 	else
2265 #endif
2266 		b = d2b(rv, &p2, &bbits);
2267 #ifdef Avoid_Underflow
2268 	p2 -= bc->scale;
2269 #endif
2270 	/* floor(log2(rv)) == bbits - 1 + p2 */
2271 	/* Check for denormal case. */
2272 	i = P - bbits;
2273 	if (i > (j = P - Emin - 1 + p2)) {
2274 #ifdef Sudden_Underflow
2275 		Bfree(b);
2276 		b = i2b(1);
2277 		p2 = Emin;
2278 		i = P - 1;
2279 #ifdef Avoid_Underflow
2280 		word0(rv) = (1 + bc->scale) << Exp_shift;
2281 #else
2282 		word0(rv) = Exp_msk1;
2283 #endif
2284 		word1(rv) = 0;
2285 #else
2286 		i = j;
2287 #endif
2288 		}
2289 #ifdef Honor_FLT_ROUNDS
2290 	if (bc->rounding != 1) {
2291 		if (i > 0)
2292 			b = lshift(b, i);
2293 		if (dsign)
2294 			b = increment(b);
2295 		}
2296 	else
2297 #endif
2298 		{
2299 		b = lshift(b, ++i);
2300 		b->x[0] |= 1;
2301 		}
2302 #ifndef Sudden_Underflow
2303  have_i:
2304 #endif
2305 	p2 -= p5 + i;
2306 	d = i2b(1);
2307 	/* Arrange for convenient computation of quotients:
2308 	 * shift left if necessary so divisor has 4 leading 0 bits.
2309 	 */
2310 	if (p5 > 0)
2311 		d = pow5mult(d, p5);
2312 	else if (p5 < 0)
2313 		b = pow5mult(b, -p5);
2314 	if (p2 > 0) {
2315 		b2 = p2;
2316 		d2 = 0;
2317 		}
2318 	else {
2319 		b2 = 0;
2320 		d2 = -p2;
2321 		}
2322 	i = dshift(d, d2);
2323 	if ((b2 += i) > 0)
2324 		b = lshift(b, b2);
2325 	if ((d2 += i) > 0)
2326 		d = lshift(d, d2);
2327 
2328 	/* Now b/d = exactly half-way between the two floating-point values */
2329 	/* on either side of the input string.  Compute first digit of b/d. */
2330 
2331 	if (!(dig = quorem(b,d))) {
2332 		b = multadd(b, 10, 0);	/* very unlikely */
2333 		dig = quorem(b,d);
2334 		}
2335 
2336 	/* Compare b/d with s0 */
2337 
2338 	for(i = 0; i < nd0; ) {
2339 		if ((dd = s0[i++] - '0' - dig))
2340 			goto ret;
2341 		if (!b->x[0] && b->wds == 1) {
2342 			if (i < nd)
2343 				dd = 1;
2344 			goto ret;
2345 			}
2346 		b = multadd(b, 10, 0);
2347 		dig = quorem(b,d);
2348 		}
2349 	for(j = bc->dp1; i++ < nd;) {
2350 		if ((dd = s0[j++] - '0' - dig))
2351 			goto ret;
2352 		if (!b->x[0] && b->wds == 1) {
2353 			if (i < nd)
2354 				dd = 1;
2355 			goto ret;
2356 			}
2357 		b = multadd(b, 10, 0);
2358 		dig = quorem(b,d);
2359 		}
2360 	if (b->x[0] || b->wds > 1)
2361 		dd = -1;
2362  ret:
2363 	Bfree(b);
2364 	Bfree(d);
2365 #ifdef Honor_FLT_ROUNDS
2366 	if (bc->rounding != 1) {
2367 		if (dd < 0) {
2368 			if (bc->rounding == 0) {
2369 				if (!dsign)
2370 					goto retlow1;
2371 				}
2372 			else if (dsign)
2373 				goto rethi1;
2374 			}
2375 		else if (dd > 0) {
2376 			if (bc->rounding == 0) {
2377 				if (dsign)
2378 					goto rethi1;
2379 				goto ret1;
2380 				}
2381 			if (!dsign)
2382 				goto rethi1;
2383 			dval(rv) += 2.*ulp(rv);
2384 			}
2385 		else {
2386 			bc->inexact = 0;
2387 			if (dsign)
2388 				goto rethi1;
2389 			}
2390 		}
2391 	else
2392 #endif
2393 	if (speccase) {
2394 		if (dd <= 0)
2395 			rv->d = 0.;
2396 		}
2397 	else if (dd < 0) {
2398 		if (!dsign)	/* does not happen for round-near */
2399 retlow1:
2400 			dval(rv) -= ulp(rv);
2401 		}
2402 	else if (dd > 0) {
2403 		if (dsign) {
2404  rethi1:
2405 			dval(rv) += ulp(rv);
2406 			}
2407 		}
2408 	else {
2409 		/* Exact half-way case:  apply round-even rule. */
2410 		if (word1(rv) & 1) {
2411 			if (dsign)
2412 				goto rethi1;
2413 			goto retlow1;
2414 			}
2415 		}
2416 
2417 #ifdef Honor_FLT_ROUNDS
2418  ret1:
2419 #endif
2420 	return;
2421 	}
2422 #endif /* NO_STRTOD_BIGCOMP */
2423 
2424  double
2425 strtod
2426 #ifdef KR_headers
2427 	(s00, se) CONST char *s00; char **se;
2428 #else
2429 	(CONST char *s00, char **se)
2430 #endif
2431 {
2432 	int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, e, e1;
2433 	int esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
2434 	CONST char *s, *s0, *s1;
2435 	double aadj, aadj1;
2436 	Long L;
2437 	U aadj2, adj, rv, rv0;
2438 	ULong y, z;
2439 	BCinfo bc;
2440 	Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
2441 #ifdef SET_INEXACT
2442 	int oldinexact;
2443 #endif
2444 #ifdef Honor_FLT_ROUNDS /*{*/
2445 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
2446 	bc.rounding = Flt_Rounds;
2447 #else /*}{*/
2448 	bc.rounding = 1;
2449 	switch(fegetround()) {
2450 	  case FE_TOWARDZERO:	bc.rounding = 0; break;
2451 	  case FE_UPWARD:	bc.rounding = 2; break;
2452 	  case FE_DOWNWARD:	bc.rounding = 3;
2453 	  }
2454 #endif /*}}*/
2455 #endif /*}*/
2456 #ifdef USE_LOCALE
2457 	CONST char *s2;
2458 #endif
2459 
2460 	sign = nz0 = nz = bc.dplen = bc.uflchk = 0;
2461 	dval(&rv) = 0.;
2462 	for(s = s00;;s++) switch(*s) {
2463 		case '-':
2464 			sign = 1;
2465 			/* no break */
2466 		case '+':
2467 			if (*++s)
2468 				goto break2;
2469 			/* no break */
2470 		case 0:
2471 			goto ret0;
2472 		case '\t':
2473 		case '\n':
2474 		case '\v':
2475 		case '\f':
2476 		case '\r':
2477 		case ' ':
2478 			continue;
2479 		default:
2480 			goto break2;
2481 		}
2482  break2:
2483 	if (*s == '0') {
2484 #ifndef NO_HEX_FP /*{*/
2485 		switch(s[1]) {
2486 		  case 'x':
2487 		  case 'X':
2488 #ifdef Honor_FLT_ROUNDS
2489 			gethex(&s, &rv, bc.rounding, sign);
2490 #else
2491 			gethex(&s, &rv, 1, sign);
2492 #endif
2493 			goto ret;
2494 		  }
2495 #endif /*}*/
2496 		nz0 = 1;
2497 		while(*++s == '0') ;
2498 		if (!*s)
2499 			goto ret;
2500 		}
2501 	s0 = s;
2502 	y = z = 0;
2503 	for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
2504 		if (nd < 9)
2505 			y = 10*y + c - '0';
2506 		else if (nd < 16)
2507 			z = 10*z + c - '0';
2508 	nd0 = nd;
2509 	bc.dp0 = bc.dp1 = s - s0;
2510 #ifdef USE_LOCALE
2511 	s1 = localeconv()->decimal_point;
2512 	if (c == *s1) {
2513 		c = '.';
2514 		if (*++s1) {
2515 			s2 = s;
2516 			for(;;) {
2517 				if (*++s2 != *s1) {
2518 					c = 0;
2519 					break;
2520 					}
2521 				if (!*++s1) {
2522 					s = s2;
2523 					break;
2524 					}
2525 				}
2526 			}
2527 		}
2528 #endif
2529 	if (c == '.') {
2530 		c = *++s;
2531 		bc.dp1 = s - s0;
2532 		bc.dplen = bc.dp1 - bc.dp0;
2533 		if (!nd) {
2534 			for(; c == '0'; c = *++s)
2535 				nz++;
2536 			if (c > '0' && c <= '9') {
2537 				s0 = s;
2538 				nf += nz;
2539 				nz = 0;
2540 				goto have_dig;
2541 				}
2542 			goto dig_done;
2543 			}
2544 		for(; c >= '0' && c <= '9'; c = *++s) {
2545  have_dig:
2546 			nz++;
2547 			if (c -= '0') {
2548 				nf += nz;
2549 				for(i = 1; i < nz; i++)
2550 					if (nd++ < 9)
2551 						y *= 10;
2552 					else if (nd <= DBL_DIG + 1)
2553 						z *= 10;
2554 				if (nd++ < 9)
2555 					y = 10*y + c;
2556 				else if (nd <= DBL_DIG + 1)
2557 					z = 10*z + c;
2558 				nz = 0;
2559 				}
2560 			}
2561 		}
2562  dig_done:
2563 	e = 0;
2564 	if (c == 'e' || c == 'E') {
2565 		if (!nd && !nz && !nz0) {
2566 			goto ret0;
2567 			}
2568 		s00 = s;
2569 		esign = 0;
2570 		switch(c = *++s) {
2571 			case '-':
2572 				esign = 1;
2573 			case '+':
2574 				c = *++s;
2575 			}
2576 		if (c >= '0' && c <= '9') {
2577 			while(c == '0')
2578 				c = *++s;
2579 			if (c > '0' && c <= '9') {
2580 				L = c - '0';
2581 				s1 = s;
2582 				while((c = *++s) >= '0' && c <= '9')
2583 					L = 10*L + c - '0';
2584 				if (s - s1 > 8 || L > 19999)
2585 					/* Avoid confusion from exponents
2586 					 * so large that e might overflow.
2587 					 */
2588 					e = 19999; /* safe for 16 bit ints */
2589 				else
2590 					e = (int)L;
2591 				if (esign)
2592 					e = -e;
2593 				}
2594 			else
2595 				e = 0;
2596 			}
2597 		else
2598 			s = s00;
2599 		}
2600 	if (!nd) {
2601 		if (!nz && !nz0) {
2602 #ifdef INFNAN_CHECK
2603 			/* Check for Nan and Infinity */
2604 			if (!bc.dplen)
2605 			 switch(c) {
2606 			  case 'i':
2607 			  case 'I':
2608 				if (match(&s,"nf")) {
2609 					--s;
2610 					if (!match(&s,"inity"))
2611 						++s;
2612 					word0(&rv) = 0x7ff00000;
2613 					word1(&rv) = 0;
2614 					goto ret;
2615 					}
2616 				break;
2617 			  case 'n':
2618 			  case 'N':
2619 				if (match(&s, "an")) {
2620 					word0(&rv) = NAN_WORD0;
2621 					word1(&rv) = NAN_WORD1;
2622 #ifndef No_Hex_NaN
2623 					if (*s == '(') /*)*/
2624 						hexnan(&rv, &s);
2625 #endif
2626 					goto ret;
2627 					}
2628 			  }
2629 #endif /* INFNAN_CHECK */
2630  ret0:
2631 			s = s00;
2632 			sign = 0;
2633 			}
2634 		goto ret;
2635 		}
2636 	bc.e0 = e1 = e -= nf;
2637 
2638 	/* Now we have nd0 digits, starting at s0, followed by a
2639 	 * decimal point, followed by nd-nd0 digits.  The number we're
2640 	 * after is the integer represented by those digits times
2641 	 * 10**e */
2642 
2643 	if (!nd0)
2644 		nd0 = nd;
2645 	k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
2646 	dval(&rv) = y;
2647 	if (k > 9) {
2648 #ifdef SET_INEXACT
2649 		if (k > DBL_DIG)
2650 			oldinexact = get_inexact();
2651 #endif
2652 		dval(&rv) = tens[k - 9] * dval(&rv) + z;
2653 		}
2654 	bd0 = 0;
2655 	if (nd <= DBL_DIG
2656 #ifndef RND_PRODQUOT
2657 #ifndef Honor_FLT_ROUNDS
2658 		&& Flt_Rounds == 1
2659 #endif
2660 #endif
2661 			) {
2662 		if (!e)
2663 			goto ret;
2664 		if (e > 0) {
2665 			if (e <= Ten_pmax) {
2666 #ifdef VAX
2667 				goto vax_ovfl_check;
2668 #else
2669 #ifdef Honor_FLT_ROUNDS
2670 				/* round correctly FLT_ROUNDS = 2 or 3 */
2671 				if (sign) {
2672 					rv.d = -rv.d;
2673 					sign = 0;
2674 					}
2675 #endif
2676 				/* rv = */ rounded_product(dval(&rv), tens[e]);
2677 				goto ret;
2678 #endif
2679 				}
2680 			i = DBL_DIG - nd;
2681 			if (e <= Ten_pmax + i) {
2682 				/* A fancier test would sometimes let us do
2683 				 * this for larger i values.
2684 				 */
2685 #ifdef Honor_FLT_ROUNDS
2686 				/* round correctly FLT_ROUNDS = 2 or 3 */
2687 				if (sign) {
2688 					rv.d = -rv.d;
2689 					sign = 0;
2690 					}
2691 #endif
2692 				e -= i;
2693 				dval(&rv) *= tens[i];
2694 #ifdef VAX
2695 				/* VAX exponent range is so narrow we must
2696 				 * worry about overflow here...
2697 				 */
2698  vax_ovfl_check:
2699 				word0(&rv) -= P*Exp_msk1;
2700 				/* rv = */ rounded_product(dval(&rv), tens[e]);
2701 				if ((word0(&rv) & Exp_mask)
2702 				 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
2703 					goto ovfl;
2704 				word0(&rv) += P*Exp_msk1;
2705 #else
2706 				/* rv = */ rounded_product(dval(&rv), tens[e]);
2707 #endif
2708 				goto ret;
2709 				}
2710 			}
2711 #ifndef Inaccurate_Divide
2712 		else if (e >= -Ten_pmax) {
2713 #ifdef Honor_FLT_ROUNDS
2714 			/* round correctly FLT_ROUNDS = 2 or 3 */
2715 			if (sign) {
2716 				rv.d = -rv.d;
2717 				sign = 0;
2718 				}
2719 #endif
2720 			/* rv = */ rounded_quotient(dval(&rv), tens[-e]);
2721 			goto ret;
2722 			}
2723 #endif
2724 		}
2725 	e1 += nd - k;
2726 
2727 #ifdef IEEE_Arith
2728 #ifdef SET_INEXACT
2729 	bc.inexact = 1;
2730 	if (k <= DBL_DIG)
2731 		oldinexact = get_inexact();
2732 #endif
2733 #ifdef Avoid_Underflow
2734 	bc.scale = 0;
2735 #endif
2736 #ifdef Honor_FLT_ROUNDS
2737 	if (bc.rounding >= 2) {
2738 		if (sign)
2739 			bc.rounding = bc.rounding == 2 ? 0 : 2;
2740 		else
2741 			if (bc.rounding != 2)
2742 				bc.rounding = 0;
2743 		}
2744 #endif
2745 #endif /*IEEE_Arith*/
2746 
2747 	/* Get starting approximation = rv * 10**e1 */
2748 
2749 	if (e1 > 0) {
2750 		if ((i = e1 & 15))
2751 			dval(&rv) *= tens[i];
2752 		if (e1 &= ~15) {
2753 			if (e1 > DBL_MAX_10_EXP) {
2754  ovfl:
2755 #ifndef NO_ERRNO
2756 				errno = ERANGE;
2757 #endif
2758 				/* Can't trust HUGE_VAL */
2759 #ifdef IEEE_Arith
2760 #ifdef Honor_FLT_ROUNDS
2761 				switch(bc.rounding) {
2762 				  case 0: /* toward 0 */
2763 				  case 3: /* toward -infinity */
2764 					word0(&rv) = Big0;
2765 					word1(&rv) = Big1;
2766 					break;
2767 				  default:
2768 					word0(&rv) = Exp_mask;
2769 					word1(&rv) = 0;
2770 				  }
2771 #else /*Honor_FLT_ROUNDS*/
2772 				word0(&rv) = Exp_mask;
2773 				word1(&rv) = 0;
2774 #endif /*Honor_FLT_ROUNDS*/
2775 #ifdef SET_INEXACT
2776 				/* set overflow bit */
2777 				dval(&rv0) = 1e300;
2778 				dval(&rv0) *= dval(&rv0);
2779 #endif
2780 #else /*IEEE_Arith*/
2781 				word0(&rv) = Big0;
2782 				word1(&rv) = Big1;
2783 #endif /*IEEE_Arith*/
2784 				goto ret;
2785 				}
2786 			e1 >>= 4;
2787 			for(j = 0; e1 > 1; j++, e1 >>= 1)
2788 				if (e1 & 1)
2789 					dval(&rv) *= bigtens[j];
2790 		/* The last multiplication could overflow. */
2791 			word0(&rv) -= P*Exp_msk1;
2792 			dval(&rv) *= bigtens[j];
2793 			if ((z = word0(&rv) & Exp_mask)
2794 			 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
2795 				goto ovfl;
2796 			if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
2797 				/* set to largest number */
2798 				/* (Can't trust DBL_MAX) */
2799 				word0(&rv) = Big0;
2800 				word1(&rv) = Big1;
2801 				}
2802 			else
2803 				word0(&rv) += P*Exp_msk1;
2804 			}
2805 		}
2806 	else if (e1 < 0) {
2807 		e1 = -e1;
2808 		if ((i = e1 & 15))
2809 			dval(&rv) /= tens[i];
2810 		if (e1 >>= 4) {
2811 			if (e1 >= 1 << n_bigtens)
2812 				goto undfl;
2813 #ifdef Avoid_Underflow
2814 			if (e1 & Scale_Bit)
2815 				bc.scale = 2*P;
2816 			for(j = 0; e1 > 0; j++, e1 >>= 1)
2817 				if (e1 & 1)
2818 					dval(&rv) *= tinytens[j];
2819 			if (bc.scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask)
2820 						>> Exp_shift)) > 0) {
2821 				/* scaled rv is denormal; clear j low bits */
2822 				if (j >= 32) {
2823 					word1(&rv) = 0;
2824 					if (j >= 53)
2825 					 word0(&rv) = (P+2)*Exp_msk1;
2826 					else
2827 					 word0(&rv) &= 0xffffffff << (j-32);
2828 					}
2829 				else
2830 					word1(&rv) &= 0xffffffff << j;
2831 				}
2832 #else
2833 			for(j = 0; e1 > 1; j++, e1 >>= 1)
2834 				if (e1 & 1)
2835 					dval(&rv) *= tinytens[j];
2836 			/* The last multiplication could underflow. */
2837 			dval(&rv0) = dval(&rv);
2838 			dval(&rv) *= tinytens[j];
2839 			if (!dval(&rv)) {
2840 				dval(&rv) = 2.*dval(&rv0);
2841 				dval(&rv) *= tinytens[j];
2842 #endif
2843 				if (!dval(&rv)) {
2844  undfl:
2845 					dval(&rv) = 0.;
2846 #ifndef NO_ERRNO
2847 					errno = ERANGE;
2848 #endif
2849 					goto ret;
2850 					}
2851 #ifndef Avoid_Underflow
2852 				word0(&rv) = Tiny0;
2853 				word1(&rv) = Tiny1;
2854 				/* The refinement below will clean
2855 				 * this approximation up.
2856 				 */
2857 				}
2858 #endif
2859 			}
2860 		}
2861 
2862 	/* Now the hard part -- adjusting rv to the correct value.*/
2863 
2864 	/* Put digits into bd: true value = bd * 10^e */
2865 
2866 	bc.nd = nd;
2867 #ifndef NO_STRTOD_BIGCOMP
2868 	bc.nd0 = nd0;	/* Only needed if nd > strtod_diglim, but done here */
2869 			/* to silence an erroneous warning about bc.nd0 */
2870 			/* possibly not being initialized. */
2871 	if (nd > strtod_diglim) {
2872 		/* ASSERT(strtod_diglim >= 18); 18 == one more than the */
2873 		/* minimum number of decimal digits to distinguish double values */
2874 		/* in IEEE arithmetic. */
2875 		i = j = 18;
2876 		if (i > nd0)
2877 			j += bc.dplen;
2878 		for(;;) {
2879 			if (--j <= bc.dp1 && j >= bc.dp0)
2880 				j = bc.dp0 - 1;
2881 			if (s0[j] != '0')
2882 				break;
2883 			--i;
2884 			}
2885 		e += nd - i;
2886 		nd = i;
2887 		if (nd0 > nd)
2888 			nd0 = nd;
2889 		if (nd < 9) { /* must recompute y */
2890 			y = 0;
2891 			for(i = 0; i < nd0; ++i)
2892 				y = 10*y + s0[i] - '0';
2893 			for(j = bc.dp1; i < nd; ++i)
2894 				y = 10*y + s0[j++] - '0';
2895 			}
2896 		}
2897 #endif
2898 	bd0 = s2b(s0, nd0, nd, y, bc.dplen);
2899 
2900 	for(;;) {
2901 		bd = Balloc(bd0->k);
2902 		Bcopy(bd, bd0);
2903 		bb = d2b(&rv, &bbe, &bbbits);	/* rv = bb * 2^bbe */
2904 		bs = i2b(1);
2905 
2906 		if (e >= 0) {
2907 			bb2 = bb5 = 0;
2908 			bd2 = bd5 = e;
2909 			}
2910 		else {
2911 			bb2 = bb5 = -e;
2912 			bd2 = bd5 = 0;
2913 			}
2914 		if (bbe >= 0)
2915 			bb2 += bbe;
2916 		else
2917 			bd2 -= bbe;
2918 		bs2 = bb2;
2919 #ifdef Honor_FLT_ROUNDS
2920 		if (bc.rounding != 1)
2921 			bs2++;
2922 #endif
2923 #ifdef Avoid_Underflow
2924 		j = bbe - bc.scale;
2925 		i = j + bbbits - 1;	/* logb(rv) */
2926 		if (i < Emin)	/* denormal */
2927 			j += P - Emin;
2928 		else
2929 			j = P + 1 - bbbits;
2930 #else /*Avoid_Underflow*/
2931 #ifdef Sudden_Underflow
2932 #ifdef IBM
2933 		j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
2934 #else
2935 		j = P + 1 - bbbits;
2936 #endif
2937 #else /*Sudden_Underflow*/
2938 		j = bbe;
2939 		i = j + bbbits - 1;	/* logb(rv) */
2940 		if (i < Emin)	/* denormal */
2941 			j += P - Emin;
2942 		else
2943 			j = P + 1 - bbbits;
2944 #endif /*Sudden_Underflow*/
2945 #endif /*Avoid_Underflow*/
2946 		bb2 += j;
2947 		bd2 += j;
2948 #ifdef Avoid_Underflow
2949 		bd2 += bc.scale;
2950 #endif
2951 		i = bb2 < bd2 ? bb2 : bd2;
2952 		if (i > bs2)
2953 			i = bs2;
2954 		if (i > 0) {
2955 			bb2 -= i;
2956 			bd2 -= i;
2957 			bs2 -= i;
2958 			}
2959 		if (bb5 > 0) {
2960 			bs = pow5mult(bs, bb5);
2961 			bb1 = mult(bs, bb);
2962 			Bfree(bb);
2963 			bb = bb1;
2964 			}
2965 		if (bb2 > 0)
2966 			bb = lshift(bb, bb2);
2967 		if (bd5 > 0)
2968 			bd = pow5mult(bd, bd5);
2969 		if (bd2 > 0)
2970 			bd = lshift(bd, bd2);
2971 		if (bs2 > 0)
2972 			bs = lshift(bs, bs2);
2973 		delta = diff(bb, bd);
2974 		bc.dsign = delta->sign;
2975 		delta->sign = 0;
2976 		i = cmp(delta, bs);
2977 #ifndef NO_STRTOD_BIGCOMP
2978 		if (bc.nd > nd && i <= 0) {
2979 			if (bc.dsign)
2980 				break;	/* Must use bigcomp(). */
2981 #ifdef Honor_FLT_ROUNDS
2982 			if (bc.rounding != 1) {
2983 				if (i < 0)
2984 					break;
2985 				}
2986 			else
2987 #endif
2988 				{
2989 				bc.nd = nd;
2990 				i = -1;	/* Discarded digits make delta smaller. */
2991 				}
2992 			}
2993 #endif
2994 #ifdef Honor_FLT_ROUNDS
2995 		if (bc.rounding != 1) {
2996 			if (i < 0) {
2997 				/* Error is less than an ulp */
2998 				if (!delta->x[0] && delta->wds <= 1) {
2999 					/* exact */
3000 #ifdef SET_INEXACT
3001 					bc.inexact = 0;
3002 #endif
3003 					break;
3004 					}
3005 				if (bc.rounding) {
3006 					if (bc.dsign) {
3007 						adj.d = 1.;
3008 						goto apply_adj;
3009 						}
3010 					}
3011 				else if (!bc.dsign) {
3012 					adj.d = -1.;
3013 					if (!word1(&rv)
3014 					 && !(word0(&rv) & Frac_mask)) {
3015 						y = word0(&rv) & Exp_mask;
3016 #ifdef Avoid_Underflow
3017 						if (!bc.scale || y > 2*P*Exp_msk1)
3018 #else
3019 						if (y)
3020 #endif
3021 						  {
3022 						  delta = lshift(delta,Log2P);
3023 						  if (cmp(delta, bs) <= 0)
3024 							adj.d = -0.5;
3025 						  }
3026 						}
3027  apply_adj:
3028 #ifdef Avoid_Underflow
3029 					if (bc.scale && (y = word0(&rv) & Exp_mask)
3030 						<= 2*P*Exp_msk1)
3031 					  word0(&adj) += (2*P+1)*Exp_msk1 - y;
3032 #else
3033 #ifdef Sudden_Underflow
3034 					if ((word0(&rv) & Exp_mask) <=
3035 							P*Exp_msk1) {
3036 						word0(&rv) += P*Exp_msk1;
3037 						dval(&rv) += adj.d*ulp(dval(&rv));
3038 						word0(&rv) -= P*Exp_msk1;
3039 						}
3040 					else
3041 #endif /*Sudden_Underflow*/
3042 #endif /*Avoid_Underflow*/
3043 					dval(&rv) += adj.d*ulp(&rv);
3044 					}
3045 				break;
3046 				}
3047 			adj.d = ratio(delta, bs);
3048 			if (adj.d < 1.)
3049 				adj.d = 1.;
3050 			if (adj.d <= 0x7ffffffe) {
3051 				/* adj = rounding ? ceil(adj) : floor(adj); */
3052 				y = adj.d;
3053 				if (y != adj.d) {
3054 					if (!((bc.rounding>>1) ^ bc.dsign))
3055 						y++;
3056 					adj.d = y;
3057 					}
3058 				}
3059 #ifdef Avoid_Underflow
3060 			if (bc.scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
3061 				word0(&adj) += (2*P+1)*Exp_msk1 - y;
3062 #else
3063 #ifdef Sudden_Underflow
3064 			if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
3065 				word0(&rv) += P*Exp_msk1;
3066 				adj.d *= ulp(dval(&rv));
3067 				if (bc.dsign)
3068 					dval(&rv) += adj.d;
3069 				else
3070 					dval(&rv) -= adj.d;
3071 				word0(&rv) -= P*Exp_msk1;
3072 				goto cont;
3073 				}
3074 #endif /*Sudden_Underflow*/
3075 #endif /*Avoid_Underflow*/
3076 			adj.d *= ulp(&rv);
3077 			if (bc.dsign) {
3078 				if (word0(&rv) == Big0 && word1(&rv) == Big1)
3079 					goto ovfl;
3080 				dval(&rv) += adj.d;
3081 				}
3082 			else
3083 				dval(&rv) -= adj.d;
3084 			goto cont;
3085 			}
3086 #endif /*Honor_FLT_ROUNDS*/
3087 
3088 		if (i < 0) {
3089 			/* Error is less than half an ulp -- check for
3090 			 * special case of mantissa a power of two.
3091 			 */
3092 			if (bc.dsign || word1(&rv) || word0(&rv) & Bndry_mask
3093 #ifdef IEEE_Arith
3094 #ifdef Avoid_Underflow
3095 			 || (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1
3096 #else
3097 			 || (word0(&rv) & Exp_mask) <= Exp_msk1
3098 #endif
3099 #endif
3100 				) {
3101 #ifdef SET_INEXACT
3102 				if (!delta->x[0] && delta->wds <= 1)
3103 					bc.inexact = 0;
3104 #endif
3105 				break;
3106 				}
3107 			if (!delta->x[0] && delta->wds <= 1) {
3108 				/* exact result */
3109 #ifdef SET_INEXACT
3110 				bc.inexact = 0;
3111 #endif
3112 				break;
3113 				}
3114 			delta = lshift(delta,Log2P);
3115 			if (cmp(delta, bs) > 0)
3116 				goto drop_down;
3117 			break;
3118 			}
3119 		if (i == 0) {
3120 			/* exactly half-way between */
3121 			if (bc.dsign) {
3122 				if ((word0(&rv) & Bndry_mask1) == Bndry_mask1
3123 				 &&  word1(&rv) == (
3124 #ifdef Avoid_Underflow
3125 			(bc.scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
3126 		? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
3127 #endif
3128 						   0xffffffff)) {
3129 					/*boundary case -- increment exponent*/
3130 					word0(&rv) = (word0(&rv) & Exp_mask)
3131 						+ Exp_msk1
3132 #ifdef IBM
3133 						| Exp_msk1 >> 4
3134 #endif
3135 						;
3136 					word1(&rv) = 0;
3137 #ifdef Avoid_Underflow
3138 					bc.dsign = 0;
3139 #endif
3140 					break;
3141 					}
3142 				}
3143 			else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) {
3144  drop_down:
3145 				/* boundary case -- decrement exponent */
3146 #ifdef Sudden_Underflow /*{{*/
3147 				L = word0(&rv) & Exp_mask;
3148 #ifdef IBM
3149 				if (L <  Exp_msk1)
3150 #else
3151 #ifdef Avoid_Underflow
3152 				if (L <= (bc.scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
3153 #else
3154 				if (L <= Exp_msk1)
3155 #endif /*Avoid_Underflow*/
3156 #endif /*IBM*/
3157 					{
3158 					if (bc.nd >nd) {
3159 						bc.uflchk = 1;
3160 						break;
3161 						}
3162 					goto undfl;
3163 					}
3164 				L -= Exp_msk1;
3165 #else /*Sudden_Underflow}{*/
3166 #ifdef Avoid_Underflow
3167 				if (bc.scale) {
3168 					L = word0(&rv) & Exp_mask;
3169 					if (L <= (2*P+1)*Exp_msk1) {
3170 						if (L > (P+2)*Exp_msk1)
3171 							/* round even ==> */
3172 							/* accept rv */
3173 							break;
3174 						/* rv = smallest denormal */
3175 						if (bc.nd >nd) {
3176 							bc.uflchk = 1;
3177 							break;
3178 							}
3179 						goto undfl;
3180 						}
3181 					}
3182 #endif /*Avoid_Underflow*/
3183 				L = (word0(&rv) & Exp_mask) - Exp_msk1;
3184 #endif /*Sudden_Underflow}}*/
3185 				word0(&rv) = L | Bndry_mask1;
3186 				word1(&rv) = 0xffffffff;
3187 #ifdef IBM
3188 				goto cont;
3189 #else
3190 				break;
3191 #endif
3192 				}
3193 #ifndef ROUND_BIASED
3194 			if (!(word1(&rv) & LSB))
3195 				break;
3196 #endif
3197 			if (bc.dsign)
3198 				dval(&rv) += ulp(&rv);
3199 #ifndef ROUND_BIASED
3200 			else {
3201 				dval(&rv) -= ulp(&rv);
3202 #ifndef Sudden_Underflow
3203 				if (!dval(&rv)) {
3204 					if (bc.nd >nd) {
3205 						bc.uflchk = 1;
3206 						break;
3207 						}
3208 					goto undfl;
3209 					}
3210 #endif
3211 				}
3212 #ifdef Avoid_Underflow
3213 			bc.dsign = 1 - bc.dsign;
3214 #endif
3215 #endif
3216 			break;
3217 			}
3218 		if ((aadj = ratio(delta, bs)) <= 2.) {
3219 			if (bc.dsign)
3220 				aadj = aadj1 = 1.;
3221 			else if (word1(&rv) || word0(&rv) & Bndry_mask) {
3222 #ifndef Sudden_Underflow
3223 				if (word1(&rv) == Tiny1 && !word0(&rv)) {
3224 					if (bc.nd >nd) {
3225 						bc.uflchk = 1;
3226 						break;
3227 						}
3228 					goto undfl;
3229 					}
3230 #endif
3231 				aadj = 1.;
3232 				aadj1 = -1.;
3233 				}
3234 			else {
3235 				/* special case -- power of FLT_RADIX to be */
3236 				/* rounded down... */
3237 
3238 				if (aadj < 2./FLT_RADIX)
3239 					aadj = 1./FLT_RADIX;
3240 				else
3241 					aadj *= 0.5;
3242 				aadj1 = -aadj;
3243 				}
3244 			}
3245 		else {
3246 			aadj *= 0.5;
3247 			aadj1 = bc.dsign ? aadj : -aadj;
3248 #ifdef Check_FLT_ROUNDS
3249 			switch(bc.rounding) {
3250 				case 2: /* towards +infinity */
3251 					aadj1 -= 0.5;
3252 					break;
3253 				case 0: /* towards 0 */
3254 				case 3: /* towards -infinity */
3255 					aadj1 += 0.5;
3256 				}
3257 #else
3258 			if (Flt_Rounds == 0)
3259 				aadj1 += 0.5;
3260 #endif /*Check_FLT_ROUNDS*/
3261 			}
3262 		y = word0(&rv) & Exp_mask;
3263 
3264 		/* Check for overflow */
3265 
3266 		if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
3267 			dval(&rv0) = dval(&rv);
3268 			word0(&rv) -= P*Exp_msk1;
3269 			adj.d = aadj1 * ulp(&rv);
3270 			dval(&rv) += adj.d;
3271 			if ((word0(&rv) & Exp_mask) >=
3272 					Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
3273 				if (word0(&rv0) == Big0 && word1(&rv0) == Big1)
3274 					goto ovfl;
3275 				word0(&rv) = Big0;
3276 				word1(&rv) = Big1;
3277 				goto cont;
3278 				}
3279 			else
3280 				word0(&rv) += P*Exp_msk1;
3281 			}
3282 		else {
3283 #ifdef Avoid_Underflow
3284 			if (bc.scale && y <= 2*P*Exp_msk1) {
3285 				if (aadj <= 0x7fffffff) {
3286 					if ((z = aadj) <= 0)
3287 						z = 1;
3288 					aadj = z;
3289 					aadj1 = bc.dsign ? aadj : -aadj;
3290 					}
3291 				dval(&aadj2) = aadj1;
3292 				word0(&aadj2) += (2*P+1)*Exp_msk1 - y;
3293 				aadj1 = dval(&aadj2);
3294 				}
3295 			adj.d = aadj1 * ulp(&rv);
3296 			dval(&rv) += adj.d;
3297 #else
3298 #ifdef Sudden_Underflow
3299 			if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
3300 				dval(&rv0) = dval(&rv);
3301 				word0(&rv) += P*Exp_msk1;
3302 				adj.d = aadj1 * ulp(&rv);
3303 				dval(&rv) += adj.d;
3304 #ifdef IBM
3305 				if ((word0(&rv) & Exp_mask) <  P*Exp_msk1)
3306 #else
3307 				if ((word0(&rv) & Exp_mask) <= P*Exp_msk1)
3308 #endif
3309 					{
3310 					if (word0(&rv0) == Tiny0
3311 					 && word1(&rv0) == Tiny1) {
3312 						if (bc.nd >nd) {
3313 							bc.uflchk = 1;
3314 							break;
3315 							}
3316 						goto undfl;
3317 						}
3318 					word0(&rv) = Tiny0;
3319 					word1(&rv) = Tiny1;
3320 					goto cont;
3321 					}
3322 				else
3323 					word0(&rv) -= P*Exp_msk1;
3324 				}
3325 			else {
3326 				adj.d = aadj1 * ulp(&rv);
3327 				dval(&rv) += adj.d;
3328 				}
3329 #else /*Sudden_Underflow*/
3330 			/* Compute adj so that the IEEE rounding rules will
3331 			 * correctly round rv + adj in some half-way cases.
3332 			 * If rv * ulp(rv) is denormalized (i.e.,
3333 			 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
3334 			 * trouble from bits lost to denormalization;
3335 			 * example: 1.2e-307 .
3336 			 */
3337 			if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
3338 				aadj1 = (double)(int)(aadj + 0.5);
3339 				if (!bc.dsign)
3340 					aadj1 = -aadj1;
3341 				}
3342 			adj.d = aadj1 * ulp(&rv);
3343 			dval(&rv) += adj.d;
3344 #endif /*Sudden_Underflow*/
3345 #endif /*Avoid_Underflow*/
3346 			}
3347 		z = word0(&rv) & Exp_mask;
3348 #ifndef SET_INEXACT
3349 		if (bc.nd == nd) {
3350 #ifdef Avoid_Underflow
3351 		if (!bc.scale)
3352 #endif
3353 		if (y == z) {
3354 			/* Can we stop now? */
3355 			L = (Long)aadj;
3356 			aadj -= L;
3357 			/* The tolerances below are conservative. */
3358 			if (bc.dsign || word1(&rv) || word0(&rv) & Bndry_mask) {
3359 				if (aadj < .4999999 || aadj > .5000001)
3360 					break;
3361 				}
3362 			else if (aadj < .4999999/FLT_RADIX)
3363 				break;
3364 			}
3365 		}
3366 #endif
3367  cont:
3368 		Bfree(bb);
3369 		Bfree(bd);
3370 		Bfree(bs);
3371 		Bfree(delta);
3372 		}
3373 	Bfree(bb);
3374 	Bfree(bd);
3375 	Bfree(bs);
3376 	Bfree(bd0);
3377 	Bfree(delta);
3378 #ifndef NO_STRTOD_BIGCOMP
3379 	if (bc.nd > nd)
3380 		bigcomp(&rv, s0, &bc);
3381 #endif
3382 #ifdef SET_INEXACT
3383 	if (bc.inexact) {
3384 		if (!oldinexact) {
3385 			word0(&rv0) = Exp_1 + (70 << Exp_shift);
3386 			word1(&rv0) = 0;
3387 			dval(&rv0) += 1.;
3388 			}
3389 		}
3390 	else if (!oldinexact)
3391 		clear_inexact();
3392 #endif
3393 #ifdef Avoid_Underflow
3394 	if (bc.scale) {
3395 		word0(&rv0) = Exp_1 - 2*P*Exp_msk1;
3396 		word1(&rv0) = 0;
3397 		dval(&rv) *= dval(&rv0);
3398 #ifndef NO_ERRNO
3399 		/* try to avoid the bug of testing an 8087 register value */
3400 #ifdef IEEE_Arith
3401 		if (!(word0(&rv) & Exp_mask))
3402 #else
3403 		if (word0(&rv) == 0 && word1(&rv) == 0)
3404 #endif
3405 			errno = ERANGE;
3406 #endif
3407 		}
3408 #endif /* Avoid_Underflow */
3409 #ifdef SET_INEXACT
3410 	if (bc.inexact && !(word0(&rv) & Exp_mask)) {
3411 		/* set underflow bit */
3412 		dval(&rv0) = 1e-300;
3413 		dval(&rv0) *= dval(&rv0);
3414 		}
3415 #endif
3416  ret:
3417 	if (se)
3418 		*se = (char *)s;
3419 	return sign ? -dval(&rv) : dval(&rv);
3420 	}
3421 
3422 #ifndef MULTIPLE_THREADS
3423  static char *dtoa_result;
3424 #endif
3425 
3426  static char *
3427 #ifdef KR_headers
3428 rv_alloc(i) int i;
3429 #else
3430 rv_alloc(int i)
3431 #endif
3432 {
3433 	int j, k, *r;
3434 
3435 	j = sizeof(ULong);
3436 	for(k = 0;
3437 		sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= (size_t)i;
3438 		j <<= 1)
3439 			k++;
3440 	r = (int*)Balloc(k);
3441 	*r = k;
3442 	return
3443 #ifndef MULTIPLE_THREADS
3444 	dtoa_result =
3445 #endif
3446 		(char *)(r+1);
3447 	}
3448 
3449  static char *
3450 #ifdef KR_headers
3451 nrv_alloc(s, rve, n) char *s, **rve; int n;
3452 #else
3453 nrv_alloc(CONST char *s, char **rve, int n)
3454 #endif
3455 {
3456 	char *rv, *t;
3457 
3458 	t = rv = rv_alloc(n);
3459 	while((*t = *s++)) t++;
3460 	if (rve)
3461 		*rve = t;
3462 	return rv;
3463 	}
3464 
3465 /* freedtoa(s) must be used to free values s returned by dtoa
3466  * when MULTIPLE_THREADS is #defined.  It should be used in all cases,
3467  * but for consistency with earlier versions of dtoa, it is optional
3468  * when MULTIPLE_THREADS is not defined.
3469  */
3470 
3471  void
3472 #ifdef KR_headers
3473 freedtoa(s) char *s;
3474 #else
3475 freedtoa(char *s)
3476 #endif
3477 {
3478 	Bigint *b = (Bigint *)((int *)s - 1);
3479 	b->maxwds = 1 << (b->k = *(int*)b);
3480 	Bfree(b);
3481 #ifndef MULTIPLE_THREADS
3482 	if (s == dtoa_result)
3483 		dtoa_result = 0;
3484 #endif
3485 	}
3486 
3487 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
3488  *
3489  * Inspired by "How to Print Floating-Point Numbers Accurately" by
3490  * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
3491  *
3492  * Modifications:
3493  *	1. Rather than iterating, we use a simple numeric overestimate
3494  *	   to determine k = floor(log10(d)).  We scale relevant
3495  *	   quantities using O(log2(k)) rather than O(k) multiplications.
3496  *	2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
3497  *	   try to generate digits strictly left to right.  Instead, we
3498  *	   compute with fewer bits and propagate the carry if necessary
3499  *	   when rounding the final digit up.  This is often faster.
3500  *	3. Under the assumption that input will be rounded nearest,
3501  *	   mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
3502  *	   That is, we allow equality in stopping tests when the
3503  *	   round-nearest rule will give the same floating-point value
3504  *	   as would satisfaction of the stopping test with strict
3505  *	   inequality.
3506  *	4. We remove common factors of powers of 2 from relevant
3507  *	   quantities.
3508  *	5. When converting floating-point integers less than 1e16,
3509  *	   we use floating-point arithmetic rather than resorting
3510  *	   to multiple-precision integers.
3511  *	6. When asked to produce fewer than 15 digits, we first try
3512  *	   to get by with floating-point arithmetic; we resort to
3513  *	   multiple-precision integer arithmetic only if we cannot
3514  *	   guarantee that the floating-point calculation has given
3515  *	   the correctly rounded result.  For k requested digits and
3516  *	   "uniformly" distributed input, the probability is
3517  *	   something like 10^(k-15) that we must resort to the Long
3518  *	   calculation.
3519  */
3520 
3521  char *
3522 dtoa
3523 #ifdef KR_headers
3524 	(dd, mode, ndigits, decpt, sign, rve)
3525 	double dd; int mode, ndigits, *decpt, *sign; char **rve;
3526 #else
3527 	(double dd, int mode, int ndigits, int *decpt, int *sign, char **rve)
3528 #endif
3529 {
3530  /*	Arguments ndigits, decpt, sign are similar to those
3531 	of ecvt and fcvt; trailing zeros are suppressed from
3532 	the returned string.  If not null, *rve is set to point
3533 	to the end of the return value.  If d is +-Infinity or NaN,
3534 	then *decpt is set to 9999.
3535 
3536 	mode:
3537 		0 ==> shortest string that yields d when read in
3538 			and rounded to nearest.
3539 		1 ==> like 0, but with Steele & White stopping rule;
3540 			e.g. with IEEE P754 arithmetic , mode 0 gives
3541 			1e23 whereas mode 1 gives 9.999999999999999e22.
3542 		2 ==> max(1,ndigits) significant digits.  This gives a
3543 			return value similar to that of ecvt, except
3544 			that trailing zeros are suppressed.
3545 		3 ==> through ndigits past the decimal point.  This
3546 			gives a return value similar to that from fcvt,
3547 			except that trailing zeros are suppressed, and
3548 			ndigits can be negative.
3549 		4,5 ==> similar to 2 and 3, respectively, but (in
3550 			round-nearest mode) with the tests of mode 0 to
3551 			possibly return a shorter string that rounds to d.
3552 			With IEEE arithmetic and compilation with
3553 			-DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
3554 			as modes 2 and 3 when FLT_ROUNDS != 1.
3555 		6-9 ==> Debugging modes similar to mode - 4:  don't try
3556 			fast floating-point estimate (if applicable).
3557 
3558 		Values of mode other than 0-9 are treated as mode 0.
3559 
3560 		Sufficient space is allocated to the return value
3561 		to hold the suppressed trailing zeros.
3562 	*/
3563 
3564 	int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
3565 		j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
3566 		spec_case, try_quick;
3567 	Long L;
3568 #ifndef Sudden_Underflow
3569 	int denorm;
3570 	ULong x;
3571 #endif
3572 	Bigint *b, *b1, *delta, *mlo, *mhi, *S;
3573 	U d2, eps, u;
3574 	double ds;
3575 	char *s, *s0;
3576 #ifdef SET_INEXACT
3577 	int inexact, oldinexact;
3578 #endif
3579 #ifdef Honor_FLT_ROUNDS /*{*/
3580 	int Rounding;
3581 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
3582 	Rounding = Flt_Rounds;
3583 #else /*}{*/
3584 	Rounding = 1;
3585 	switch(fegetround()) {
3586 	  case FE_TOWARDZERO:	Rounding = 0; break;
3587 	  case FE_UPWARD:	Rounding = 2; break;
3588 	  case FE_DOWNWARD:	Rounding = 3;
3589 	  }
3590 #endif /*}}*/
3591 #endif /*}*/
3592 
3593 #ifndef MULTIPLE_THREADS
3594 	if (dtoa_result) {
3595 		freedtoa(dtoa_result);
3596 		dtoa_result = 0;
3597 		}
3598 #endif
3599 
3600 	u.d = dd;
3601 	if (word0(&u) & Sign_bit) {
3602 		/* set sign for everything, including 0's and NaNs */
3603 		*sign = 1;
3604 		word0(&u) &= ~Sign_bit;	/* clear sign bit */
3605 		}
3606 	else
3607 		*sign = 0;
3608 
3609 #if defined(IEEE_Arith) + defined(VAX)
3610 #ifdef IEEE_Arith
3611 	if ((word0(&u) & Exp_mask) == Exp_mask)
3612 #else
3613 	if (word0(&u)  == 0x8000)
3614 #endif
3615 		{
3616 		/* Infinity or NaN */
3617 		*decpt = 9999;
3618 #ifdef IEEE_Arith
3619 		if (!word1(&u) && !(word0(&u) & 0xfffff))
3620 			return nrv_alloc("Infinity", rve, 8);
3621 #endif
3622 		return nrv_alloc("NaN", rve, 3);
3623 		}
3624 #endif
3625 #ifdef IBM
3626 	dval(&u) += 0; /* normalize */
3627 #endif
3628 	if (!dval(&u)) {
3629 		*decpt = 1;
3630 		return nrv_alloc("0", rve, 1);
3631 		}
3632 
3633 #ifdef SET_INEXACT
3634 	try_quick = oldinexact = get_inexact();
3635 	inexact = 1;
3636 #endif
3637 #ifdef Honor_FLT_ROUNDS
3638 	if (Rounding >= 2) {
3639 		if (*sign)
3640 			Rounding = Rounding == 2 ? 0 : 2;
3641 		else
3642 			if (Rounding != 2)
3643 				Rounding = 0;
3644 		}
3645 #endif
3646 
3647 	b = d2b(&u, &be, &bbits);
3648 #ifdef Sudden_Underflow
3649 	i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
3650 #else
3651 	if ((i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1)))) {
3652 #endif
3653 		dval(&d2) = dval(&u);
3654 		word0(&d2) &= Frac_mask1;
3655 		word0(&d2) |= Exp_11;
3656 #ifdef IBM
3657 		if (j = 11 - hi0bits(word0(&d2) & Frac_mask))
3658 			dval(&d2) /= 1 << j;
3659 #endif
3660 
3661 		/* log(x)	~=~ log(1.5) + (x-1.5)/1.5
3662 		 * log10(x)	 =  log(x) / log(10)
3663 		 *		~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
3664 		 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
3665 		 *
3666 		 * This suggests computing an approximation k to log10(d) by
3667 		 *
3668 		 * k = (i - Bias)*0.301029995663981
3669 		 *	+ ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
3670 		 *
3671 		 * We want k to be too large rather than too small.
3672 		 * The error in the first-order Taylor series approximation
3673 		 * is in our favor, so we just round up the constant enough
3674 		 * to compensate for any error in the multiplication of
3675 		 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
3676 		 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
3677 		 * adding 1e-13 to the constant term more than suffices.
3678 		 * Hence we adjust the constant term to 0.1760912590558.
3679 		 * (We could get a more accurate k by invoking log10,
3680 		 *  but this is probably not worthwhile.)
3681 		 */
3682 
3683 		i -= Bias;
3684 #ifdef IBM
3685 		i <<= 2;
3686 		i += j;
3687 #endif
3688 #ifndef Sudden_Underflow
3689 		denorm = 0;
3690 		}
3691 	else {
3692 		/* d is denormalized */
3693 
3694 		i = bbits + be + (Bias + (P-1) - 1);
3695 		x = i > 32  ? word0(&u) << (64 - i) | word1(&u) >> (i - 32)
3696 			    : word1(&u) << (32 - i);
3697 		dval(&d2) = x;
3698 		word0(&d2) -= 31*Exp_msk1; /* adjust exponent */
3699 		i -= (Bias + (P-1) - 1) + 1;
3700 		denorm = 1;
3701 		}
3702 #endif
3703 	ds = (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
3704 	k = (int)ds;
3705 	if (ds < 0. && ds != k)
3706 		k--;	/* want k = floor(ds) */
3707 	k_check = 1;
3708 	if (k >= 0 && k <= Ten_pmax) {
3709 		if (dval(&u) < tens[k])
3710 			k--;
3711 		k_check = 0;
3712 		}
3713 	j = bbits - i - 1;
3714 	if (j >= 0) {
3715 		b2 = 0;
3716 		s2 = j;
3717 		}
3718 	else {
3719 		b2 = -j;
3720 		s2 = 0;
3721 		}
3722 	if (k >= 0) {
3723 		b5 = 0;
3724 		s5 = k;
3725 		s2 += k;
3726 		}
3727 	else {
3728 		b2 -= k;
3729 		b5 = -k;
3730 		s5 = 0;
3731 		}
3732 	if (mode < 0 || mode > 9)
3733 		mode = 0;
3734 
3735 #ifndef SET_INEXACT
3736 #ifdef Check_FLT_ROUNDS
3737 	try_quick = Rounding == 1;
3738 #else
3739 	try_quick = 1;
3740 #endif
3741 #endif /*SET_INEXACT*/
3742 
3743 	if (mode > 5) {
3744 		mode -= 4;
3745 		try_quick = 0;
3746 		}
3747 	leftright = 1;
3748 	ilim = ilim1 = -1;	/* Values for cases 0 and 1; done here to */
3749 				/* silence erroneous "gcc -Wall" warning. */
3750 	switch(mode) {
3751 		case 0:
3752 		case 1:
3753 			i = 18;
3754 			ndigits = 0;
3755 			break;
3756 		case 2:
3757 			leftright = 0;
3758 			/* no break */
3759 		case 4:
3760 			if (ndigits <= 0)
3761 				ndigits = 1;
3762 			ilim = ilim1 = i = ndigits;
3763 			break;
3764 		case 3:
3765 			leftright = 0;
3766 			/* no break */
3767 		case 5:
3768 			i = ndigits + k + 1;
3769 			ilim = i;
3770 			ilim1 = i - 1;
3771 			if (i <= 0)
3772 				i = 1;
3773 		}
3774 	s = s0 = rv_alloc(i);
3775 
3776 #ifdef Honor_FLT_ROUNDS
3777 	if (mode > 1 && Rounding != 1)
3778 		leftright = 0;
3779 #endif
3780 
3781 	if (ilim >= 0 && ilim <= Quick_max && try_quick) {
3782 
3783 		/* Try to get by with floating-point arithmetic. */
3784 
3785 		i = 0;
3786 		dval(&d2) = dval(&u);
3787 		k0 = k;
3788 		ilim0 = ilim;
3789 		ieps = 2; /* conservative */
3790 		if (k > 0) {
3791 			ds = tens[k&0xf];
3792 			j = k >> 4;
3793 			if (j & Bletch) {
3794 				/* prevent overflows */
3795 				j &= Bletch - 1;
3796 				dval(&u) /= bigtens[n_bigtens-1];
3797 				ieps++;
3798 				}
3799 			for(; j; j >>= 1, i++)
3800 				if (j & 1) {
3801 					ieps++;
3802 					ds *= bigtens[i];
3803 					}
3804 			dval(&u) /= ds;
3805 			}
3806 		else if ((j1 = -k)) {
3807 			dval(&u) *= tens[j1 & 0xf];
3808 			for(j = j1 >> 4; j; j >>= 1, i++)
3809 				if (j & 1) {
3810 					ieps++;
3811 					dval(&u) *= bigtens[i];
3812 					}
3813 			}
3814 		if (k_check && dval(&u) < 1. && ilim > 0) {
3815 			if (ilim1 <= 0)
3816 				goto fast_failed;
3817 			ilim = ilim1;
3818 			k--;
3819 			dval(&u) *= 10.;
3820 			ieps++;
3821 			}
3822 		dval(&eps) = ieps*dval(&u) + 7.;
3823 		word0(&eps) -= (P-1)*Exp_msk1;
3824 		if (ilim == 0) {
3825 			S = mhi = 0;
3826 			dval(&u) -= 5.;
3827 			if (dval(&u) > dval(&eps))
3828 				goto one_digit;
3829 			if (dval(&u) < -dval(&eps))
3830 				goto no_digits;
3831 			goto fast_failed;
3832 			}
3833 #ifndef No_leftright
3834 		if (leftright) {
3835 			/* Use Steele & White method of only
3836 			 * generating digits needed.
3837 			 */
3838 			dval(&eps) = 0.5/tens[ilim-1] - dval(&eps);
3839 			for(i = 0;;) {
3840 				L = dval(&u);
3841 				dval(&u) -= L;
3842 				*s++ = '0' + (int)L;
3843 				if (dval(&u) < dval(&eps))
3844 					goto ret1;
3845 				if (1. - dval(&u) < dval(&eps))
3846 					goto bump_up;
3847 				if (++i >= ilim)
3848 					break;
3849 				dval(&eps) *= 10.;
3850 				dval(&u) *= 10.;
3851 				}
3852 			}
3853 		else {
3854 #endif
3855 			/* Generate ilim digits, then fix them up. */
3856 			dval(&eps) *= tens[ilim-1];
3857 			for(i = 1;; i++, dval(&u) *= 10.) {
3858 				L = (Long)(dval(&u));
3859 				if (!(dval(&u) -= L))
3860 					ilim = i;
3861 				*s++ = '0' + (int)L;
3862 				if (i == ilim) {
3863 					if (dval(&u) > 0.5 + dval(&eps))
3864 						goto bump_up;
3865 					else if (dval(&u) < 0.5 - dval(&eps)) {
3866 						while(*--s == '0') {}
3867 						s++;
3868 						goto ret1;
3869 						}
3870 					break;
3871 					}
3872 				}
3873 #ifndef No_leftright
3874 			}
3875 #endif
3876  fast_failed:
3877 		s = s0;
3878 		dval(&u) = dval(&d2);
3879 		k = k0;
3880 		ilim = ilim0;
3881 		}
3882 
3883 	/* Do we have a "small" integer? */
3884 
3885 	if (be >= 0 && k <= Int_max) {
3886 		/* Yes. */
3887 		ds = tens[k];
3888 		if (ndigits < 0 && ilim <= 0) {
3889 			S = mhi = 0;
3890 			if (ilim < 0 || dval(&u) <= 5*ds)
3891 				goto no_digits;
3892 			goto one_digit;
3893 			}
3894 		for(i = 1;; i++, dval(&u) *= 10.) {
3895 			L = (Long)(dval(&u) / ds);
3896 			dval(&u) -= L*ds;
3897 #ifdef Check_FLT_ROUNDS
3898 			/* If FLT_ROUNDS == 2, L will usually be high by 1 */
3899 			if (dval(&u) < 0) {
3900 				L--;
3901 				dval(&u) += ds;
3902 				}
3903 #endif
3904 			*s++ = '0' + (int)L;
3905 			if (!dval(&u)) {
3906 #ifdef SET_INEXACT
3907 				inexact = 0;
3908 #endif
3909 				break;
3910 				}
3911 			if (i == ilim) {
3912 #ifdef Honor_FLT_ROUNDS
3913 				if (mode > 1)
3914 				switch(Rounding) {
3915 				  case 0: goto ret1;
3916 				  case 2: goto bump_up;
3917 				  }
3918 #endif
3919 				dval(&u) += dval(&u);
3920 				if (dval(&u) > ds || (dval(&u) == ds && L & 1)) {
3921  bump_up:
3922 					while(*--s == '9')
3923 						if (s == s0) {
3924 							k++;
3925 							*s = '0';
3926 							break;
3927 							}
3928 					++*s++;
3929 					}
3930 				break;
3931 				}
3932 			}
3933 		goto ret1;
3934 		}
3935 
3936 	m2 = b2;
3937 	m5 = b5;
3938 	mhi = mlo = 0;
3939 	if (leftright) {
3940 		i =
3941 #ifndef Sudden_Underflow
3942 			denorm ? be + (Bias + (P-1) - 1 + 1) :
3943 #endif
3944 #ifdef IBM
3945 			1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
3946 #else
3947 			1 + P - bbits;
3948 #endif
3949 		b2 += i;
3950 		s2 += i;
3951 		mhi = i2b(1);
3952 		}
3953 	if (m2 > 0 && s2 > 0) {
3954 		i = m2 < s2 ? m2 : s2;
3955 		b2 -= i;
3956 		m2 -= i;
3957 		s2 -= i;
3958 		}
3959 	if (b5 > 0) {
3960 		if (leftright) {
3961 			if (m5 > 0) {
3962 				mhi = pow5mult(mhi, m5);
3963 				b1 = mult(mhi, b);
3964 				Bfree(b);
3965 				b = b1;
3966 				}
3967 			if ((j = b5 - m5))
3968 				b = pow5mult(b, j);
3969 			}
3970 		else
3971 			b = pow5mult(b, b5);
3972 		}
3973 	S = i2b(1);
3974 	if (s5 > 0)
3975 		S = pow5mult(S, s5);
3976 
3977 	/* Check for special case that d is a normalized power of 2. */
3978 
3979 	spec_case = 0;
3980 	if ((mode < 2 || leftright)
3981 #ifdef Honor_FLT_ROUNDS
3982 			&& Rounding == 1
3983 #endif
3984 				) {
3985 		if (!word1(&u) && !(word0(&u) & Bndry_mask)
3986 #ifndef Sudden_Underflow
3987 		 && word0(&u) & (Exp_mask & ~Exp_msk1)
3988 #endif
3989 				) {
3990 			/* The special case */
3991 			b2 += Log2P;
3992 			s2 += Log2P;
3993 			spec_case = 1;
3994 			}
3995 		}
3996 
3997 	/* Arrange for convenient computation of quotients:
3998 	 * shift left if necessary so divisor has 4 leading 0 bits.
3999 	 *
4000 	 * Perhaps we should just compute leading 28 bits of S once
4001 	 * and for all and pass them and a shift to quorem, so it
4002 	 * can do shifts and ors to compute the numerator for q.
4003 	 */
4004 #ifdef Pack_32
4005 	if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f))
4006 		i = 32 - i;
4007 #define iInc 28
4008 #else
4009 	if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf)
4010 		i = 16 - i;
4011 #define iInc 12
4012 #endif
4013 	i = dshift(S, s2);
4014 	b2 += i;
4015 	m2 += i;
4016 	s2 += i;
4017 	if (b2 > 0)
4018 		b = lshift(b, b2);
4019 	if (s2 > 0)
4020 		S = lshift(S, s2);
4021 	if (k_check) {
4022 		if (cmp(b,S) < 0) {
4023 			k--;
4024 			b = multadd(b, 10, 0);	/* we botched the k estimate */
4025 			if (leftright)
4026 				mhi = multadd(mhi, 10, 0);
4027 			ilim = ilim1;
4028 			}
4029 		}
4030 	if (ilim <= 0 && (mode == 3 || mode == 5)) {
4031 		if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
4032 			/* no digits, fcvt style */
4033  no_digits:
4034 			k = -1 - ndigits;
4035 			goto ret;
4036 			}
4037  one_digit:
4038 		*s++ = '1';
4039 		k++;
4040 		goto ret;
4041 		}
4042 	if (leftright) {
4043 		if (m2 > 0)
4044 			mhi = lshift(mhi, m2);
4045 
4046 		/* Compute mlo -- check for special case
4047 		 * that d is a normalized power of 2.
4048 		 */
4049 
4050 		mlo = mhi;
4051 		if (spec_case) {
4052 			mhi = Balloc(mhi->k);
4053 			Bcopy(mhi, mlo);
4054 			mhi = lshift(mhi, Log2P);
4055 			}
4056 
4057 		for(i = 1;;i++) {
4058 			dig = quorem(b,S) + '0';
4059 			/* Do we yet have the shortest decimal string
4060 			 * that will round to d?
4061 			 */
4062 			j = cmp(b, mlo);
4063 			delta = diff(S, mhi);
4064 			j1 = delta->sign ? 1 : cmp(b, delta);
4065 			Bfree(delta);
4066 #ifndef ROUND_BIASED
4067 			if (j1 == 0 && mode != 1 && !(word1(&u) & 1)
4068 #ifdef Honor_FLT_ROUNDS
4069 				&& Rounding >= 1
4070 #endif
4071 								   ) {
4072 				if (dig == '9')
4073 					goto round_9_up;
4074 				if (j > 0)
4075 					dig++;
4076 #ifdef SET_INEXACT
4077 				else if (!b->x[0] && b->wds <= 1)
4078 					inexact = 0;
4079 #endif
4080 				*s++ = dig;
4081 				goto ret;
4082 				}
4083 #endif
4084 			if (j < 0 || (j == 0 && mode != 1
4085 #ifndef ROUND_BIASED
4086 							&& !(word1(&u) & 1)
4087 #endif
4088 					)) {
4089 				if (!b->x[0] && b->wds <= 1) {
4090 #ifdef SET_INEXACT
4091 					inexact = 0;
4092 #endif
4093 					goto accept_dig;
4094 					}
4095 #ifdef Honor_FLT_ROUNDS
4096 				if (mode > 1)
4097 				 switch(Rounding) {
4098 				  case 0: goto accept_dig;
4099 				  case 2: goto keep_dig;
4100 				  }
4101 #endif /*Honor_FLT_ROUNDS*/
4102 				if (j1 > 0) {
4103 					b = lshift(b, 1);
4104 					j1 = cmp(b, S);
4105 					if ((j1 > 0 || (j1 == 0 && dig & 1))
4106 					&& dig++ == '9')
4107 						goto round_9_up;
4108 					}
4109  accept_dig:
4110 				*s++ = dig;
4111 				goto ret;
4112 				}
4113 			if (j1 > 0) {
4114 #ifdef Honor_FLT_ROUNDS
4115 				if (!Rounding)
4116 					goto accept_dig;
4117 #endif
4118 				if (dig == '9') { /* possible if i == 1 */
4119  round_9_up:
4120 					*s++ = '9';
4121 					goto roundoff;
4122 					}
4123 				*s++ = dig + 1;
4124 				goto ret;
4125 				}
4126 #ifdef Honor_FLT_ROUNDS
4127  keep_dig:
4128 #endif
4129 			*s++ = dig;
4130 			if (i == ilim)
4131 				break;
4132 			b = multadd(b, 10, 0);
4133 			if (mlo == mhi)
4134 				mlo = mhi = multadd(mhi, 10, 0);
4135 			else {
4136 				mlo = multadd(mlo, 10, 0);
4137 				mhi = multadd(mhi, 10, 0);
4138 				}
4139 			}
4140 		}
4141 	else
4142 		for(i = 1;; i++) {
4143 			*s++ = dig = quorem(b,S) + '0';
4144 			if (!b->x[0] && b->wds <= 1) {
4145 #ifdef SET_INEXACT
4146 				inexact = 0;
4147 #endif
4148 				goto ret;
4149 				}
4150 			if (i >= ilim)
4151 				break;
4152 			b = multadd(b, 10, 0);
4153 			}
4154 
4155 	/* Round off last digit */
4156 
4157 #ifdef Honor_FLT_ROUNDS
4158 	switch(Rounding) {
4159 	  case 0: goto trimzeros;
4160 	  case 2: goto roundoff;
4161 	  }
4162 #endif
4163 	b = lshift(b, 1);
4164 	j = cmp(b, S);
4165 	if (j > 0 || (j == 0 && dig & 1)) {
4166  roundoff:
4167 		while(*--s == '9')
4168 			if (s == s0) {
4169 				k++;
4170 				*s++ = '1';
4171 				goto ret;
4172 				}
4173 		++*s++;
4174 		}
4175 	else {
4176 #ifdef Honor_FLT_ROUNDS
4177  trimzeros:
4178 #endif
4179 		while(*--s == '0') {}
4180 		s++;
4181 		}
4182  ret:
4183 	Bfree(S);
4184 	if (mhi) {
4185 		if (mlo && mlo != mhi)
4186 			Bfree(mlo);
4187 		Bfree(mhi);
4188 		}
4189  ret1:
4190 #ifdef SET_INEXACT
4191 	if (inexact) {
4192 		if (!oldinexact) {
4193 			word0(&u) = Exp_1 + (70 << Exp_shift);
4194 			word1(&u) = 0;
4195 			dval(&u) += 1.;
4196 			}
4197 		}
4198 	else if (!oldinexact)
4199 		clear_inexact();
4200 #endif
4201 	Bfree(b);
4202 	*s = 0;
4203 	*decpt = k + 1;
4204 	if (rve)
4205 		*rve = s;
4206 	return s0;
4207 	}
4208 
4209 }  // namespace dmg_fp
4210