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