xref: /php-src/Zend/zend_strtod.c (revision d28f6e69)
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 #ifndef CONST
296 #ifdef KR_headers
297 #define CONST /* blank */
298 #else
299 #define CONST const
300 #endif
301 #endif
302 
303 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
304 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
305 #endif
306 
307 typedef union { double d; ULong L[2]; } U;
308 
309 #ifdef IEEE_8087
310 #define word0(x) (x)->L[1]
311 #define word1(x) (x)->L[0]
312 #else
313 #define word0(x) (x)->L[0]
314 #define word1(x) (x)->L[1]
315 #endif
316 #define dval(x) (x)->d
317 
318 #ifndef STRTOD_DIGLIM
319 #define STRTOD_DIGLIM 40
320 #endif
321 
322 #ifdef DIGLIM_DEBUG
323 extern int strtod_diglim;
324 #else
325 #define strtod_diglim STRTOD_DIGLIM
326 #endif
327 
328 /* The following definition of Storeinc is appropriate for MIPS processors.
329  * An alternative that might be better on some machines is
330  * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
331  */
332 #if defined(IEEE_8087) + defined(VAX) + defined(__arm__)
333 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
334 ((unsigned short *)a)[0] = (unsigned short)c, a++)
335 #else
336 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
337 ((unsigned short *)a)[1] = (unsigned short)c, a++)
338 #endif
339 
340 /* #define P DBL_MANT_DIG */
341 /* Ten_pmax = floor(P*log(2)/log(5)) */
342 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
343 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
344 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
345 
346 #ifdef IEEE_Arith
347 #define Exp_shift  20
348 #define Exp_shift1 20
349 #define Exp_msk1    0x100000
350 #define Exp_msk11   0x100000
351 #define Exp_mask  0x7ff00000
352 #define P 53
353 #define Nbits 53
354 #define Bias 1023
355 #define Emax 1023
356 #define Emin (-1022)
357 #define Exp_1  0x3ff00000
358 #define Exp_11 0x3ff00000
359 #define Ebits 11
360 #define Frac_mask  0xfffff
361 #define Frac_mask1 0xfffff
362 #define Ten_pmax 22
363 #define Bletch 0x10
364 #define Bndry_mask  0xfffff
365 #define Bndry_mask1 0xfffff
366 #define LSB 1
367 #define Sign_bit 0x80000000
368 #define Log2P 1
369 #define Tiny0 0
370 #define Tiny1 1
371 #define Quick_max 14
372 #define Int_max 14
373 #ifndef NO_IEEE_Scale
374 #define Avoid_Underflow
375 #ifdef Flush_Denorm	/* debugging option */
376 #undef Sudden_Underflow
377 #endif
378 #endif
379 
380 #ifndef Flt_Rounds
381 #ifdef FLT_ROUNDS
382 #define Flt_Rounds FLT_ROUNDS
383 #else
384 #define Flt_Rounds 1
385 #endif
386 #endif /*Flt_Rounds*/
387 
388 #ifdef Honor_FLT_ROUNDS
389 #undef Check_FLT_ROUNDS
390 #define Check_FLT_ROUNDS
391 #else
392 #define Rounding Flt_Rounds
393 #endif
394 
395 #else /* ifndef IEEE_Arith */
396 #undef Check_FLT_ROUNDS
397 #undef Honor_FLT_ROUNDS
398 #undef SET_INEXACT
399 #undef  Sudden_Underflow
400 #define Sudden_Underflow
401 #ifdef IBM
402 #undef Flt_Rounds
403 #define Flt_Rounds 0
404 #define Exp_shift  24
405 #define Exp_shift1 24
406 #define Exp_msk1   0x1000000
407 #define Exp_msk11  0x1000000
408 #define Exp_mask  0x7f000000
409 #define P 14
410 #define Nbits 56
411 #define Bias 65
412 #define Emax 248
413 #define Emin (-260)
414 #define Exp_1  0x41000000
415 #define Exp_11 0x41000000
416 #define Ebits 8	/* exponent has 7 bits, but 8 is the right value in b2d */
417 #define Frac_mask  0xffffff
418 #define Frac_mask1 0xffffff
419 #define Bletch 4
420 #define Ten_pmax 22
421 #define Bndry_mask  0xefffff
422 #define Bndry_mask1 0xffffff
423 #define LSB 1
424 #define Sign_bit 0x80000000
425 #define Log2P 4
426 #define Tiny0 0x100000
427 #define Tiny1 0
428 #define Quick_max 14
429 #define Int_max 15
430 #else /* VAX */
431 #undef Flt_Rounds
432 #define Flt_Rounds 1
433 #define Exp_shift  23
434 #define Exp_shift1 7
435 #define Exp_msk1    0x80
436 #define Exp_msk11   0x800000
437 #define Exp_mask  0x7f80
438 #define P 56
439 #define Nbits 56
440 #define Bias 129
441 #define Emax 126
442 #define Emin (-129)
443 #define Exp_1  0x40800000
444 #define Exp_11 0x4080
445 #define Ebits 8
446 #define Frac_mask  0x7fffff
447 #define Frac_mask1 0xffff007f
448 #define Ten_pmax 24
449 #define Bletch 2
450 #define Bndry_mask  0xffff007f
451 #define Bndry_mask1 0xffff007f
452 #define LSB 0x10000
453 #define Sign_bit 0x8000
454 #define Log2P 1
455 #define Tiny0 0x80
456 #define Tiny1 0
457 #define Quick_max 15
458 #define Int_max 15
459 #endif /* IBM, VAX */
460 #endif /* IEEE_Arith */
461 
462 #ifndef IEEE_Arith
463 #define ROUND_BIASED
464 #else
465 #ifdef ROUND_BIASED_without_Round_Up
466 #undef  ROUND_BIASED
467 #define ROUND_BIASED
468 #endif
469 #endif
470 
471 #ifdef RND_PRODQUOT
472 #define rounded_product(a,b) a = rnd_prod(a, b)
473 #define rounded_quotient(a,b) a = rnd_quot(a, b)
474 #ifdef KR_headers
475 extern double rnd_prod(), rnd_quot();
476 #else
477 extern double rnd_prod(double, double), rnd_quot(double, double);
478 #endif
479 #else
480 #define rounded_product(a,b) a *= b
481 #define rounded_quotient(a,b) a /= b
482 #endif
483 
484 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
485 #define Big1 0xffffffff
486 
487 #ifndef Pack_32
488 #define Pack_32
489 #endif
490 
491 typedef struct BCinfo BCinfo;
492  struct
493 BCinfo { int dp0, dp1, dplen, dsign, e0, inexact, nd, nd0, rounding, scale, uflchk; };
494 
495 #ifdef KR_headers
496 #define FFFFFFFF ((((unsigned long)0xffff)<<16)|(unsigned long)0xffff)
497 #else
498 #define FFFFFFFF 0xffffffffUL
499 #endif
500 
501 #ifdef NO_LONG_LONG
502 #undef ULLong
503 #ifdef Just_16
504 #undef Pack_32
505 /* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
506  * This makes some inner loops simpler and sometimes saves work
507  * during multiplications, but it often seems to make things slightly
508  * slower.  Hence the default is now to store 32 bits per Long.
509  */
510 #endif
511 #else	/* long long available */
512 #ifndef Llong
513 #define Llong long long
514 #endif
515 #ifndef ULLong
516 #define ULLong unsigned Llong
517 #endif
518 #endif /* NO_LONG_LONG */
519 
520 #ifndef MULTIPLE_THREADS
521 #define ACQUIRE_DTOA_LOCK(n)	/*nothing*/
522 #define FREE_DTOA_LOCK(n)	/*nothing*/
523 #endif
524 
525 #define Kmax 7
526 
527  struct
528 Bigint {
529 	struct Bigint *next;
530 	int k, maxwds, sign, wds;
531 	ULong x[1];
532 	};
533 
534  typedef struct Bigint Bigint;
535 
536  static Bigint *freelist[Kmax+1];
537 
538 static void destroy_freelist(void);
539 static void free_p5s(void);
540 
541 #ifdef ZTS
542 static MUTEX_T dtoa_mutex;
543 static MUTEX_T pow5mult_mutex;
544 #endif /* ZTS */
545 
zend_startup_strtod(void)546 ZEND_API int zend_startup_strtod(void) /* {{{ */
547 {
548 #ifdef ZTS
549 	dtoa_mutex = tsrm_mutex_alloc();
550 	pow5mult_mutex = tsrm_mutex_alloc();
551 #endif
552 	return 1;
553 }
554 /* }}} */
zend_shutdown_strtod(void)555 ZEND_API int zend_shutdown_strtod(void) /* {{{ */
556 {
557 	destroy_freelist();
558 	free_p5s();
559 
560 #ifdef ZTS
561 	tsrm_mutex_free(dtoa_mutex);
562 	dtoa_mutex = NULL;
563 
564 	tsrm_mutex_free(pow5mult_mutex);
565 	pow5mult_mutex = NULL;
566 #endif
567 	return 1;
568 }
569 /* }}} */
570 
571  static Bigint *
Balloc(k)572 Balloc
573 #ifdef KR_headers
574 	(k) int k;
575 #else
576 	(int k)
577 #endif
578 {
579 	int x;
580 	Bigint *rv;
581 #ifndef Omit_Private_Memory
582 	unsigned int len;
583 #endif
584 
585 	ACQUIRE_DTOA_LOCK(0);
586 	/* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */
587 	/* but this case seems very unlikely. */
588 	if (k <= Kmax && (rv = freelist[k]))
589 		freelist[k] = rv->next;
590 	else {
591 		x = 1 << k;
592 #ifdef Omit_Private_Memory
593 		rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong));
594 		if (!rv) {
595 			FREE_DTOA_LOCK(0);
596 			zend_error_noreturn(E_ERROR, "Balloc() failed to allocate memory");
597 		}
598 #else
599 		len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
600 			/sizeof(double);
601 		if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) {
602 			rv = (Bigint*)pmem_next;
603 			pmem_next += len;
604 			}
605 		else
606 			rv = (Bigint*)MALLOC(len*sizeof(double));
607 			if (!rv) {
608 				FREE_DTOA_LOCK(0);
609 				zend_error_noreturn(E_ERROR, "Balloc() failed to allocate memory");
610 			}
611 #endif
612 		rv->k = k;
613 		rv->maxwds = x;
614 		}
615 	FREE_DTOA_LOCK(0);
616 	rv->sign = rv->wds = 0;
617 	return rv;
618 	}
619 
620  static void
Bfree(v)621 Bfree
622 #ifdef KR_headers
623 	(v) Bigint *v;
624 #else
625 	(Bigint *v)
626 #endif
627 {
628 	if (v) {
629 		if (v->k > Kmax)
630 #ifdef FREE
631 			FREE((void*)v);
632 #else
633 			free((void*)v);
634 #endif
635 		else {
636 			ACQUIRE_DTOA_LOCK(0);
637 			v->next = freelist[v->k];
638 			freelist[v->k] = v;
639 			FREE_DTOA_LOCK(0);
640 			}
641 		}
642 	}
643 
644 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
645 y->wds*sizeof(Long) + 2*sizeof(int))
646 
647  static Bigint *
multadd(b,m,a)648 multadd
649 #ifdef KR_headers
650 	(b, m, a) Bigint *b; int m, a;
651 #else
652 	(Bigint *b, int m, int a)	/* multiply by m and add a */
653 #endif
654 {
655 	int i, wds;
656 #ifdef ULLong
657 	ULong *x;
658 	ULLong carry, y;
659 #else
660 	ULong carry, *x, y;
661 #ifdef Pack_32
662 	ULong xi, z;
663 #endif
664 #endif
665 	Bigint *b1;
666 
667 	wds = b->wds;
668 	x = b->x;
669 	i = 0;
670 	carry = a;
671 	do {
672 #ifdef ULLong
673 		y = *x * (ULLong)m + carry;
674 		carry = y >> 32;
675 		*x++ = y & FFFFFFFF;
676 #else
677 #ifdef Pack_32
678 		xi = *x;
679 		y = (xi & 0xffff) * m + carry;
680 		z = (xi >> 16) * m + (y >> 16);
681 		carry = z >> 16;
682 		*x++ = (z << 16) + (y & 0xffff);
683 #else
684 		y = *x * m + carry;
685 		carry = y >> 16;
686 		*x++ = y & 0xffff;
687 #endif
688 #endif
689 		}
690 		while(++i < wds);
691 	if (carry) {
692 		if (wds >= b->maxwds) {
693 			b1 = Balloc(b->k+1);
694 			Bcopy(b1, b);
695 			Bfree(b);
696 			b = b1;
697 			}
698 		b->x[wds++] = carry;
699 		b->wds = wds;
700 		}
701 	return b;
702 	}
703 
704  static Bigint *
s2b(s,nd0,nd,y9,dplen)705 s2b
706 #ifdef KR_headers
707 	(s, nd0, nd, y9, dplen) CONST char *s; int nd0, nd, dplen; ULong y9;
708 #else
709 	(const char *s, int nd0, int nd, ULong y9, int dplen)
710 #endif
711 {
712 	Bigint *b;
713 	int i, k;
714 	Long x, y;
715 
716 	x = (nd + 8) / 9;
717 	for(k = 0, y = 1; x > y; y <<= 1, k++) ;
718 #ifdef Pack_32
719 	b = Balloc(k);
720 	b->x[0] = y9;
721 	b->wds = 1;
722 #else
723 	b = Balloc(k+1);
724 	b->x[0] = y9 & 0xffff;
725 	b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
726 #endif
727 
728 	i = 9;
729 	if (9 < nd0) {
730 		s += 9;
731 		do b = multadd(b, 10, *s++ - '0');
732 			while(++i < nd0);
733 		s += dplen;
734 		}
735 	else
736 		s += dplen + 9;
737 	for(; i < nd; i++)
738 		b = multadd(b, 10, *s++ - '0');
739 	return b;
740 	}
741 
742  static int
hi0bits(x)743 hi0bits
744 #ifdef KR_headers
745 	(x) ULong x;
746 #else
747 	(ULong x)
748 #endif
749 {
750 	int k = 0;
751 
752 	if (!(x & 0xffff0000)) {
753 		k = 16;
754 		x <<= 16;
755 		}
756 	if (!(x & 0xff000000)) {
757 		k += 8;
758 		x <<= 8;
759 		}
760 	if (!(x & 0xf0000000)) {
761 		k += 4;
762 		x <<= 4;
763 		}
764 	if (!(x & 0xc0000000)) {
765 		k += 2;
766 		x <<= 2;
767 		}
768 	if (!(x & 0x80000000)) {
769 		k++;
770 		if (!(x & 0x40000000))
771 			return 32;
772 		}
773 	return k;
774 	}
775 
776  static int
lo0bits(y)777 lo0bits
778 #ifdef KR_headers
779 	(y) ULong *y;
780 #else
781 	(ULong *y)
782 #endif
783 {
784 	int k;
785 	ULong x = *y;
786 
787 	if (x & 7) {
788 		if (x & 1)
789 			return 0;
790 		if (x & 2) {
791 			*y = x >> 1;
792 			return 1;
793 			}
794 		*y = x >> 2;
795 		return 2;
796 		}
797 	k = 0;
798 	if (!(x & 0xffff)) {
799 		k = 16;
800 		x >>= 16;
801 		}
802 	if (!(x & 0xff)) {
803 		k += 8;
804 		x >>= 8;
805 		}
806 	if (!(x & 0xf)) {
807 		k += 4;
808 		x >>= 4;
809 		}
810 	if (!(x & 0x3)) {
811 		k += 2;
812 		x >>= 2;
813 		}
814 	if (!(x & 1)) {
815 		k++;
816 		x >>= 1;
817 		if (!x)
818 			return 32;
819 		}
820 	*y = x;
821 	return k;
822 	}
823 
824  static Bigint *
i2b(i)825 i2b
826 #ifdef KR_headers
827 	(i) int i;
828 #else
829 	(int i)
830 #endif
831 {
832 	Bigint *b;
833 
834 	b = Balloc(1);
835 	b->x[0] = i;
836 	b->wds = 1;
837 	return b;
838 	}
839 
840  static Bigint *
mult(a,b)841 mult
842 #ifdef KR_headers
843 	(a, b) Bigint *a, *b;
844 #else
845 	(Bigint *a, Bigint *b)
846 #endif
847 {
848 	Bigint *c;
849 	int k, wa, wb, wc;
850 	ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
851 	ULong y;
852 #ifdef ULLong
853 	ULLong carry, z;
854 #else
855 	ULong carry, z;
856 #ifdef Pack_32
857 	ULong z2;
858 #endif
859 #endif
860 
861 	if (a->wds < b->wds) {
862 		c = a;
863 		a = b;
864 		b = c;
865 		}
866 	k = a->k;
867 	wa = a->wds;
868 	wb = b->wds;
869 	wc = wa + wb;
870 	if (wc > a->maxwds)
871 		k++;
872 	c = Balloc(k);
873 	for(x = c->x, xa = x + wc; x < xa; x++)
874 		*x = 0;
875 	xa = a->x;
876 	xae = xa + wa;
877 	xb = b->x;
878 	xbe = xb + wb;
879 	xc0 = c->x;
880 #ifdef ULLong
881 	for(; xb < xbe; xc0++) {
882 		if ((y = *xb++)) {
883 			x = xa;
884 			xc = xc0;
885 			carry = 0;
886 			do {
887 				z = *x++ * (ULLong)y + *xc + carry;
888 				carry = z >> 32;
889 				*xc++ = z & FFFFFFFF;
890 				}
891 				while(x < xae);
892 			*xc = carry;
893 			}
894 		}
895 #else
896 #ifdef Pack_32
897 	for(; xb < xbe; xb++, xc0++) {
898 		if (y = *xb & 0xffff) {
899 			x = xa;
900 			xc = xc0;
901 			carry = 0;
902 			do {
903 				z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
904 				carry = z >> 16;
905 				z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
906 				carry = z2 >> 16;
907 				Storeinc(xc, z2, z);
908 				}
909 				while(x < xae);
910 			*xc = carry;
911 			}
912 		if (y = *xb >> 16) {
913 			x = xa;
914 			xc = xc0;
915 			carry = 0;
916 			z2 = *xc;
917 			do {
918 				z = (*x & 0xffff) * y + (*xc >> 16) + carry;
919 				carry = z >> 16;
920 				Storeinc(xc, z, z2);
921 				z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
922 				carry = z2 >> 16;
923 				}
924 				while(x < xae);
925 			*xc = z2;
926 			}
927 		}
928 #else
929 	for(; xb < xbe; xc0++) {
930 		if (y = *xb++) {
931 			x = xa;
932 			xc = xc0;
933 			carry = 0;
934 			do {
935 				z = *x++ * y + *xc + carry;
936 				carry = z >> 16;
937 				*xc++ = z & 0xffff;
938 				}
939 				while(x < xae);
940 			*xc = carry;
941 			}
942 		}
943 #endif
944 #endif
945 	for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
946 	c->wds = wc;
947 	return c;
948 	}
949 
950  static Bigint *p5s;
951 
952  static Bigint *
pow5mult(b,k)953 pow5mult
954 #ifdef KR_headers
955 	(b, k) Bigint *b; int k;
956 #else
957 	(Bigint *b, int k)
958 #endif
959 {
960 	Bigint *b1, *p5, *p51;
961 	int i;
962 	static int p05[3] = { 5, 25, 125 };
963 
964 	if ((i = k & 3))
965 		b = multadd(b, p05[i-1], 0);
966 
967 	if (!(k >>= 2))
968 		return b;
969 	if (!(p5 = p5s)) {
970 		/* first time */
971 #ifdef MULTIPLE_THREADS
972 		ACQUIRE_DTOA_LOCK(1);
973 		if (!(p5 = p5s)) {
974 			p5 = p5s = i2b(625);
975 			p5->next = 0;
976 			}
977 		FREE_DTOA_LOCK(1);
978 #else
979 		p5 = p5s = i2b(625);
980 		p5->next = 0;
981 #endif
982 		}
983 	for(;;) {
984 		if (k & 1) {
985 			b1 = mult(b, p5);
986 			Bfree(b);
987 			b = b1;
988 			}
989 		if (!(k >>= 1))
990 			break;
991 		if (!(p51 = p5->next)) {
992 #ifdef MULTIPLE_THREADS
993 			ACQUIRE_DTOA_LOCK(1);
994 			if (!(p51 = p5->next)) {
995 				p51 = p5->next = mult(p5,p5);
996 				p51->next = 0;
997 				}
998 			FREE_DTOA_LOCK(1);
999 #else
1000 			p51 = p5->next = mult(p5,p5);
1001 			p51->next = 0;
1002 #endif
1003 			}
1004 		p5 = p51;
1005 		}
1006 	return b;
1007 	}
1008 
1009  static Bigint *
lshift(b,k)1010 lshift
1011 #ifdef KR_headers
1012 	(b, k) Bigint *b; int k;
1013 #else
1014 	(Bigint *b, int k)
1015 #endif
1016 {
1017 	int i, k1, n, n1;
1018 	Bigint *b1;
1019 	ULong *x, *x1, *xe, z;
1020 
1021 #ifdef Pack_32
1022 	n = k >> 5;
1023 #else
1024 	n = k >> 4;
1025 #endif
1026 	k1 = b->k;
1027 	n1 = n + b->wds + 1;
1028 	for(i = b->maxwds; n1 > i; i <<= 1)
1029 		k1++;
1030 	b1 = Balloc(k1);
1031 	x1 = b1->x;
1032 	for(i = 0; i < n; i++)
1033 		*x1++ = 0;
1034 	x = b->x;
1035 	xe = x + b->wds;
1036 #ifdef Pack_32
1037 	if (k &= 0x1f) {
1038 		k1 = 32 - k;
1039 		z = 0;
1040 		do {
1041 			*x1++ = *x << k | z;
1042 			z = *x++ >> k1;
1043 			}
1044 			while(x < xe);
1045 		if ((*x1 = z))
1046 			++n1;
1047 		}
1048 #else
1049 	if (k &= 0xf) {
1050 		k1 = 16 - k;
1051 		z = 0;
1052 		do {
1053 			*x1++ = *x << k  & 0xffff | z;
1054 			z = *x++ >> k1;
1055 			}
1056 			while(x < xe);
1057 		if (*x1 = z)
1058 			++n1;
1059 		}
1060 #endif
1061 	else do
1062 		*x1++ = *x++;
1063 		while(x < xe);
1064 	b1->wds = n1 - 1;
1065 	Bfree(b);
1066 	return b1;
1067 	}
1068 
1069  static int
cmp(a,b)1070 cmp
1071 #ifdef KR_headers
1072 	(a, b) Bigint *a, *b;
1073 #else
1074 	(Bigint *a, Bigint *b)
1075 #endif
1076 {
1077 	ULong *xa, *xa0, *xb, *xb0;
1078 	int i, j;
1079 
1080 	i = a->wds;
1081 	j = b->wds;
1082 #ifdef DEBUG
1083 	if (i > 1 && !a->x[i-1])
1084 		Bug("cmp called with a->x[a->wds-1] == 0");
1085 	if (j > 1 && !b->x[j-1])
1086 		Bug("cmp called with b->x[b->wds-1] == 0");
1087 #endif
1088 	if (i -= j)
1089 		return i;
1090 	xa0 = a->x;
1091 	xa = xa0 + j;
1092 	xb0 = b->x;
1093 	xb = xb0 + j;
1094 	for(;;) {
1095 		if (*--xa != *--xb)
1096 			return *xa < *xb ? -1 : 1;
1097 		if (xa <= xa0)
1098 			break;
1099 		}
1100 	return 0;
1101 	}
1102 
1103  static Bigint *
diff(a,b)1104 diff
1105 #ifdef KR_headers
1106 	(a, b) Bigint *a, *b;
1107 #else
1108 	(Bigint *a, Bigint *b)
1109 #endif
1110 {
1111 	Bigint *c;
1112 	int i, wa, wb;
1113 	ULong *xa, *xae, *xb, *xbe, *xc;
1114 #ifdef ULLong
1115 	ULLong borrow, y;
1116 #else
1117 	ULong borrow, y;
1118 #ifdef Pack_32
1119 	ULong z;
1120 #endif
1121 #endif
1122 
1123 	i = cmp(a,b);
1124 	if (!i) {
1125 		c = Balloc(0);
1126 		c->wds = 1;
1127 		c->x[0] = 0;
1128 		return c;
1129 		}
1130 	if (i < 0) {
1131 		c = a;
1132 		a = b;
1133 		b = c;
1134 		i = 1;
1135 		}
1136 	else
1137 		i = 0;
1138 	c = Balloc(a->k);
1139 	c->sign = i;
1140 	wa = a->wds;
1141 	xa = a->x;
1142 	xae = xa + wa;
1143 	wb = b->wds;
1144 	xb = b->x;
1145 	xbe = xb + wb;
1146 	xc = c->x;
1147 	borrow = 0;
1148 #ifdef ULLong
1149 	do {
1150 		y = (ULLong)*xa++ - *xb++ - borrow;
1151 		borrow = y >> 32 & (ULong)1;
1152 		*xc++ = y & FFFFFFFF;
1153 		}
1154 		while(xb < xbe);
1155 	while(xa < xae) {
1156 		y = *xa++ - borrow;
1157 		borrow = y >> 32 & (ULong)1;
1158 		*xc++ = y & FFFFFFFF;
1159 		}
1160 #else
1161 #ifdef Pack_32
1162 	do {
1163 		y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
1164 		borrow = (y & 0x10000) >> 16;
1165 		z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
1166 		borrow = (z & 0x10000) >> 16;
1167 		Storeinc(xc, z, y);
1168 		}
1169 		while(xb < xbe);
1170 	while(xa < xae) {
1171 		y = (*xa & 0xffff) - borrow;
1172 		borrow = (y & 0x10000) >> 16;
1173 		z = (*xa++ >> 16) - borrow;
1174 		borrow = (z & 0x10000) >> 16;
1175 		Storeinc(xc, z, y);
1176 		}
1177 #else
1178 	do {
1179 		y = *xa++ - *xb++ - borrow;
1180 		borrow = (y & 0x10000) >> 16;
1181 		*xc++ = y & 0xffff;
1182 		}
1183 		while(xb < xbe);
1184 	while(xa < xae) {
1185 		y = *xa++ - borrow;
1186 		borrow = (y & 0x10000) >> 16;
1187 		*xc++ = y & 0xffff;
1188 		}
1189 #endif
1190 #endif
1191 	while(!*--xc)
1192 		wa--;
1193 	c->wds = wa;
1194 	return c;
1195 	}
1196 
1197  static double
ulp(x)1198 ulp
1199 #ifdef KR_headers
1200 	(x) U *x;
1201 #else
1202 	(U *x)
1203 #endif
1204 {
1205 	Long L;
1206 	U u;
1207 
1208 	L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
1209 #ifndef Avoid_Underflow
1210 #ifndef Sudden_Underflow
1211 	if (L > 0) {
1212 #endif
1213 #endif
1214 #ifdef IBM
1215 		L |= Exp_msk1 >> 4;
1216 #endif
1217 		word0(&u) = L;
1218 		word1(&u) = 0;
1219 #ifndef Avoid_Underflow
1220 #ifndef Sudden_Underflow
1221 		}
1222 	else {
1223 		L = -L >> Exp_shift;
1224 		if (L < Exp_shift) {
1225 			word0(&u) = 0x80000 >> L;
1226 			word1(&u) = 0;
1227 			}
1228 		else {
1229 			word0(&u) = 0;
1230 			L -= Exp_shift;
1231 			word1(&u) = L >= 31 ? 1 : 1 << 31 - L;
1232 			}
1233 		}
1234 #endif
1235 #endif
1236 	return dval(&u);
1237 	}
1238 
1239  static double
b2d(a,e)1240 b2d
1241 #ifdef KR_headers
1242 	(a, e) Bigint *a; int *e;
1243 #else
1244 	(Bigint *a, int *e)
1245 #endif
1246 {
1247 	ULong *xa, *xa0, w, y, z;
1248 	int k;
1249 	U d;
1250 #ifdef VAX
1251 	ULong d0, d1;
1252 #else
1253 #define d0 word0(&d)
1254 #define d1 word1(&d)
1255 #endif
1256 
1257 	xa0 = a->x;
1258 	xa = xa0 + a->wds;
1259 	y = *--xa;
1260 #ifdef DEBUG
1261 	if (!y) Bug("zero y in b2d");
1262 #endif
1263 	k = hi0bits(y);
1264 	*e = 32 - k;
1265 #ifdef Pack_32
1266 	if (k < Ebits) {
1267 		d0 = Exp_1 | y >> (Ebits - k);
1268 		w = xa > xa0 ? *--xa : 0;
1269 		d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
1270 		goto ret_d;
1271 		}
1272 	z = xa > xa0 ? *--xa : 0;
1273 	if (k -= Ebits) {
1274 		d0 = Exp_1 | y << k | z >> (32 - k);
1275 		y = xa > xa0 ? *--xa : 0;
1276 		d1 = z << k | y >> (32 - k);
1277 		}
1278 	else {
1279 		d0 = Exp_1 | y;
1280 		d1 = z;
1281 		}
1282 #else
1283 	if (k < Ebits + 16) {
1284 		z = xa > xa0 ? *--xa : 0;
1285 		d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
1286 		w = xa > xa0 ? *--xa : 0;
1287 		y = xa > xa0 ? *--xa : 0;
1288 		d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
1289 		goto ret_d;
1290 		}
1291 	z = xa > xa0 ? *--xa : 0;
1292 	w = xa > xa0 ? *--xa : 0;
1293 	k -= Ebits + 16;
1294 	d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
1295 	y = xa > xa0 ? *--xa : 0;
1296 	d1 = w << k + 16 | y << k;
1297 #endif
1298  ret_d:
1299 #ifdef VAX
1300 	word0(&d) = d0 >> 16 | d0 << 16;
1301 	word1(&d) = d1 >> 16 | d1 << 16;
1302 #else
1303 #undef d0
1304 #undef d1
1305 #endif
1306 	return dval(&d);
1307 	}
1308 
1309  static Bigint *
d2b(d,e,bits)1310 d2b
1311 #ifdef KR_headers
1312 	(d, e, bits) U *d; int *e, *bits;
1313 #else
1314 	(U *d, int *e, int *bits)
1315 #endif
1316 {
1317 	Bigint *b;
1318 	int de, k;
1319 	ULong *x, y, z;
1320 #ifndef Sudden_Underflow
1321 	int i;
1322 #endif
1323 #ifdef VAX
1324 	ULong d0, d1;
1325 	d0 = word0(d) >> 16 | word0(d) << 16;
1326 	d1 = word1(d) >> 16 | word1(d) << 16;
1327 #else
1328 #define d0 word0(d)
1329 #define d1 word1(d)
1330 #endif
1331 
1332 #ifdef Pack_32
1333 	b = Balloc(1);
1334 #else
1335 	b = Balloc(2);
1336 #endif
1337 	x = b->x;
1338 
1339 	z = d0 & Frac_mask;
1340 	d0 &= 0x7fffffff;	/* clear sign bit, which we ignore */
1341 #ifdef Sudden_Underflow
1342 	de = (int)(d0 >> Exp_shift);
1343 #ifndef IBM
1344 	z |= Exp_msk11;
1345 #endif
1346 #else
1347 	if ((de = (int)(d0 >> Exp_shift)))
1348 		z |= Exp_msk1;
1349 #endif
1350 #ifdef Pack_32
1351 	if ((y = d1)) {
1352 		if ((k = lo0bits(&y))) {
1353 			x[0] = y | z << (32 - k);
1354 			z >>= k;
1355 			}
1356 		else
1357 			x[0] = y;
1358 #ifndef Sudden_Underflow
1359 		i =
1360 #endif
1361 		    b->wds = (x[1] = z) ? 2 : 1;
1362 		}
1363 	else {
1364 		k = lo0bits(&z);
1365 		x[0] = z;
1366 #ifndef Sudden_Underflow
1367 		i =
1368 #endif
1369 		    b->wds = 1;
1370 		k += 32;
1371 		}
1372 #else
1373 	if (y = d1) {
1374 		if (k = lo0bits(&y))
1375 			if (k >= 16) {
1376 				x[0] = y | z << 32 - k & 0xffff;
1377 				x[1] = z >> k - 16 & 0xffff;
1378 				x[2] = z >> k;
1379 				i = 2;
1380 				}
1381 			else {
1382 				x[0] = y & 0xffff;
1383 				x[1] = y >> 16 | z << 16 - k & 0xffff;
1384 				x[2] = z >> k & 0xffff;
1385 				x[3] = z >> k+16;
1386 				i = 3;
1387 				}
1388 		else {
1389 			x[0] = y & 0xffff;
1390 			x[1] = y >> 16;
1391 			x[2] = z & 0xffff;
1392 			x[3] = z >> 16;
1393 			i = 3;
1394 			}
1395 		}
1396 	else {
1397 #ifdef DEBUG
1398 		if (!z)
1399 			Bug("Zero passed to d2b");
1400 #endif
1401 		k = lo0bits(&z);
1402 		if (k >= 16) {
1403 			x[0] = z;
1404 			i = 0;
1405 			}
1406 		else {
1407 			x[0] = z & 0xffff;
1408 			x[1] = z >> 16;
1409 			i = 1;
1410 			}
1411 		k += 32;
1412 		}
1413 	while(!x[i])
1414 		--i;
1415 	b->wds = i + 1;
1416 #endif
1417 #ifndef Sudden_Underflow
1418 	if (de) {
1419 #endif
1420 #ifdef IBM
1421 		*e = (de - Bias - (P-1) << 2) + k;
1422 		*bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
1423 #else
1424 		*e = de - Bias - (P-1) + k;
1425 		*bits = P - k;
1426 #endif
1427 #ifndef Sudden_Underflow
1428 		}
1429 	else {
1430 		*e = de - Bias - (P-1) + 1 + k;
1431 #ifdef Pack_32
1432 		*bits = 32*i - hi0bits(x[i-1]);
1433 #else
1434 		*bits = (i+2)*16 - hi0bits(x[i]);
1435 #endif
1436 		}
1437 #endif
1438 	return b;
1439 	}
1440 #undef d0
1441 #undef d1
1442 
1443  static double
ratio(a,b)1444 ratio
1445 #ifdef KR_headers
1446 	(a, b) Bigint *a, *b;
1447 #else
1448 	(Bigint *a, Bigint *b)
1449 #endif
1450 {
1451 	U da, db;
1452 	int k, ka, kb;
1453 
1454 	dval(&da) = b2d(a, &ka);
1455 	dval(&db) = b2d(b, &kb);
1456 #ifdef Pack_32
1457 	k = ka - kb + 32*(a->wds - b->wds);
1458 #else
1459 	k = ka - kb + 16*(a->wds - b->wds);
1460 #endif
1461 #ifdef IBM
1462 	if (k > 0) {
1463 		word0(&da) += (k >> 2)*Exp_msk1;
1464 		if (k &= 3)
1465 			dval(&da) *= 1 << k;
1466 		}
1467 	else {
1468 		k = -k;
1469 		word0(&db) += (k >> 2)*Exp_msk1;
1470 		if (k &= 3)
1471 			dval(&db) *= 1 << k;
1472 		}
1473 #else
1474 	if (k > 0)
1475 		word0(&da) += k*Exp_msk1;
1476 	else {
1477 		k = -k;
1478 		word0(&db) += k*Exp_msk1;
1479 		}
1480 #endif
1481 	return dval(&da) / dval(&db);
1482 	}
1483 
1484  static CONST double
1485 tens[] = {
1486 		1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1487 		1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1488 		1e20, 1e21, 1e22
1489 #ifdef VAX
1490 		, 1e23, 1e24
1491 #endif
1492 		};
1493 
1494  static CONST double
1495 #ifdef IEEE_Arith
1496 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1497 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
1498 #ifdef Avoid_Underflow
1499 		9007199254740992.*9007199254740992.e-256
1500 		/* = 2^106 * 1e-256 */
1501 #else
1502 		1e-256
1503 #endif
1504 		};
1505 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
1506 /* flag unnecessarily.  It leads to a song and dance at the end of strtod. */
1507 #define Scale_Bit 0x10
1508 #define n_bigtens 5
1509 #else
1510 #ifdef IBM
1511 bigtens[] = { 1e16, 1e32, 1e64 };
1512 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
1513 #define n_bigtens 3
1514 #else
1515 bigtens[] = { 1e16, 1e32 };
1516 static CONST double tinytens[] = { 1e-16, 1e-32 };
1517 #define n_bigtens 2
1518 #endif
1519 #endif
1520 
1521 #undef Need_Hexdig
1522 #ifdef INFNAN_CHECK
1523 #ifndef No_Hex_NaN
1524 #define Need_Hexdig
1525 #endif
1526 #endif
1527 
1528 #ifndef Need_Hexdig
1529 #ifndef NO_HEX_FP
1530 #define Need_Hexdig
1531 #endif
1532 #endif
1533 
1534 #ifdef Need_Hexdig /*{*/
1535 #if 0
1536 static unsigned char hexdig[256];
1537 
1538  static void
1539 htinit(unsigned char *h, unsigned char *s, int inc)
1540 {
1541 	int i, j;
1542 	for(i = 0; (j = s[i]) !=0; i++)
1543 		h[j] = i + inc;
1544 	}
1545 
1546  static void
1547 hexdig_init(void)	/* Use of hexdig_init omitted 20121220 to avoid a */
1548 			/* race condition when multiple threads are used. */
1549 {
1550 #define USC (unsigned char *)
1551 	htinit(hexdig, USC "0123456789", 0x10);
1552 	htinit(hexdig, USC "abcdef", 0x10 + 10);
1553 	htinit(hexdig, USC "ABCDEF", 0x10 + 10);
1554 	}
1555 #else
1556 static const unsigned char hexdig[256] = {
1557 	0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
1558 	0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
1559 	0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
1560 	16,17,18,19,20,21,22,23,24,25,0,0,0,0,0,0,
1561 	0,26,27,28,29,30,31,0,0,0,0,0,0,0,0,0,
1562 	0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
1563 	0,26,27,28,29,30,31,0,0,0,0,0,0,0,0,0,
1564 	0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
1565 	0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
1566 	0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
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 	0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
1571 	0,0,0,0,0,0,0,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 	};
1574 #endif
1575 #endif /* } Need_Hexdig */
1576 
1577 #ifdef INFNAN_CHECK
1578 
1579 #ifndef NAN_WORD0
1580 #define NAN_WORD0 0x7ff80000
1581 #endif
1582 
1583 #ifndef NAN_WORD1
1584 #define NAN_WORD1 0
1585 #endif
1586 
1587  static int
match(sp,t)1588 match
1589 #ifdef KR_headers
1590 	(sp, t) char **sp, *t;
1591 #else
1592 	(const char **sp, const char *t)
1593 #endif
1594 {
1595 	int c, d;
1596 	CONST char *s = *sp;
1597 
1598 	while((d = *t++)) {
1599 		if ((c = *++s) >= 'A' && c <= 'Z')
1600 			c += 'a' - 'A';
1601 		if (c != d)
1602 			return 0;
1603 		}
1604 	*sp = s + 1;
1605 	return 1;
1606 	}
1607 
1608 #ifndef No_Hex_NaN
1609  static void
hexnan(rvp,sp)1610 hexnan
1611 #ifdef KR_headers
1612 	(rvp, sp) U *rvp; CONST char **sp;
1613 #else
1614 	(U *rvp, const char **sp)
1615 #endif
1616 {
1617 	ULong c, x[2];
1618 	CONST char *s;
1619 	int c1, havedig, udx0, xshift;
1620 
1621 	/**** if (!hexdig['0']) hexdig_init(); ****/
1622 	x[0] = x[1] = 0;
1623 	havedig = xshift = 0;
1624 	udx0 = 1;
1625 	s = *sp;
1626 	/* allow optional initial 0x or 0X */
1627 	while((c = *(CONST unsigned char*)(s+1)) && c <= ' ')
1628 		++s;
1629 	if (s[1] == '0' && (s[2] == 'x' || s[2] == 'X'))
1630 		s += 2;
1631 	while((c = *(CONST unsigned char*)++s)) {
1632 		if ((c1 = hexdig[c]))
1633 			c  = c1 & 0xf;
1634 		else if (c <= ' ') {
1635 			if (udx0 && havedig) {
1636 				udx0 = 0;
1637 				xshift = 1;
1638 				}
1639 			continue;
1640 			}
1641 #ifdef GDTOA_NON_PEDANTIC_NANCHECK
1642 		else if (/*(*/ c == ')' && havedig) {
1643 			*sp = s + 1;
1644 			break;
1645 			}
1646 		else
1647 			return;	/* invalid form: don't change *sp */
1648 #else
1649 		else {
1650 			do {
1651 				if (/*(*/ c == ')') {
1652 					*sp = s + 1;
1653 					break;
1654 					}
1655 				} while((c = *++s));
1656 			break;
1657 			}
1658 #endif
1659 		havedig = 1;
1660 		if (xshift) {
1661 			xshift = 0;
1662 			x[0] = x[1];
1663 			x[1] = 0;
1664 			}
1665 		if (udx0)
1666 			x[0] = (x[0] << 4) | (x[1] >> 28);
1667 		x[1] = (x[1] << 4) | c;
1668 		}
1669 	if ((x[0] &= 0xfffff) || x[1]) {
1670 		word0(rvp) = Exp_mask | x[0];
1671 		word1(rvp) = x[1];
1672 		}
1673 	}
1674 #endif /*No_Hex_NaN*/
1675 #endif /* INFNAN_CHECK */
1676 
1677 #ifdef Pack_32
1678 #define ULbits 32
1679 #define kshift 5
1680 #define kmask 31
1681 #else
1682 #define ULbits 16
1683 #define kshift 4
1684 #define kmask 15
1685 #endif
1686 
1687 #if !defined(NO_HEX_FP) || defined(Honor_FLT_ROUNDS) /*{*/
1688  static Bigint *
1689 #ifdef KR_headers
increment(b)1690 increment(b) Bigint *b;
1691 #else
1692 increment(Bigint *b)
1693 #endif
1694 {
1695 	ULong *x, *xe;
1696 	Bigint *b1;
1697 
1698 	x = b->x;
1699 	xe = x + b->wds;
1700 	do {
1701 		if (*x < (ULong)0xffffffffL) {
1702 			++*x;
1703 			return b;
1704 			}
1705 		*x++ = 0;
1706 		} while(x < xe);
1707 	{
1708 		if (b->wds >= b->maxwds) {
1709 			b1 = Balloc(b->k+1);
1710 			Bcopy(b1,b);
1711 			Bfree(b);
1712 			b = b1;
1713 			}
1714 		b->x[b->wds++] = 1;
1715 		}
1716 	return b;
1717 	}
1718 
1719 #endif /*}*/
1720 
1721 #ifndef NO_HEX_FP /*{*/
1722 
1723  static void
1724 #ifdef KR_headers
rshift(b,k)1725 rshift(b, k) Bigint *b; int k;
1726 #else
1727 rshift(Bigint *b, int k)
1728 #endif
1729 {
1730 	ULong *x, *x1, *xe, y;
1731 	int n;
1732 
1733 	x = x1 = b->x;
1734 	n = k >> kshift;
1735 	if (n < b->wds) {
1736 		xe = x + b->wds;
1737 		x += n;
1738 		if (k &= kmask) {
1739 			n = 32 - k;
1740 			y = *x++ >> k;
1741 			while(x < xe) {
1742 				*x1++ = (y | (*x << n)) & 0xffffffff;
1743 				y = *x++ >> k;
1744 				}
1745 			if ((*x1 = y) !=0)
1746 				x1++;
1747 			}
1748 		else
1749 			while(x < xe)
1750 				*x1++ = *x++;
1751 		}
1752 	if ((b->wds = x1 - b->x) == 0)
1753 		b->x[0] = 0;
1754 	}
1755 
1756  static ULong
1757 #ifdef KR_headers
any_on(b,k)1758 any_on(b, k) Bigint *b; int k;
1759 #else
1760 any_on(Bigint *b, int k)
1761 #endif
1762 {
1763 	int n, nwds;
1764 	ULong *x, *x0, x1, x2;
1765 
1766 	x = b->x;
1767 	nwds = b->wds;
1768 	n = k >> kshift;
1769 	if (n > nwds)
1770 		n = nwds;
1771 	else if (n < nwds && (k &= kmask)) {
1772 		x1 = x2 = x[n];
1773 		x1 >>= k;
1774 		x1 <<= k;
1775 		if (x1 != x2)
1776 			return 1;
1777 		}
1778 	x0 = x;
1779 	x += n;
1780 	while(x > x0)
1781 		if (*--x)
1782 			return 1;
1783 	return 0;
1784 	}
1785 
1786 enum {	/* rounding values: same as FLT_ROUNDS */
1787 	Round_zero = 0,
1788 	Round_near = 1,
1789 	Round_up = 2,
1790 	Round_down = 3
1791 	};
1792 
1793  void
1794 #ifdef KR_headers
gethex(sp,rvp,rounding,sign)1795 gethex(sp, rvp, rounding, sign)
1796 	CONST char **sp; U *rvp; int rounding, sign;
1797 #else
1798 gethex( CONST char **sp, U *rvp, int rounding, int sign)
1799 #endif
1800 {
1801 	Bigint *b;
1802 	CONST unsigned char *decpt, *s0, *s, *s1;
1803 	Long e, e1;
1804 	ULong L, lostbits, *x;
1805 	int big, denorm, esign, havedig, k, n, nbits, up, zret;
1806 #ifdef IBM
1807 	int j;
1808 #endif
1809 	enum {
1810 #ifdef IEEE_Arith /*{{*/
1811 		emax = 0x7fe - Bias - P + 1,
1812 		emin = Emin - P + 1
1813 #else /*}{*/
1814 		emin = Emin - P,
1815 #ifdef VAX
1816 		emax = 0x7ff - Bias - P + 1
1817 #endif
1818 #ifdef IBM
1819 		emax = 0x7f - Bias - P
1820 #endif
1821 #endif /*}}*/
1822 		};
1823 #ifdef USE_LOCALE
1824 	int i;
1825 #ifdef NO_LOCALE_CACHE
1826 	const unsigned char *decimalpoint = (unsigned char*)
1827 		localeconv()->decimal_point;
1828 #else
1829 	const unsigned char *decimalpoint;
1830 	static unsigned char *decimalpoint_cache;
1831 	if (!(s0 = decimalpoint_cache)) {
1832 		s0 = (unsigned char*)localeconv()->decimal_point;
1833 		if ((decimalpoint_cache = (unsigned char*)
1834 				MALLOC(strlen((CONST char*)s0) + 1))) {
1835 			strcpy((char*)decimalpoint_cache, (CONST char*)s0);
1836 			s0 = decimalpoint_cache;
1837 			}
1838 		}
1839 	decimalpoint = s0;
1840 #endif
1841 #endif
1842 
1843 	/**** if (!hexdig['0']) hexdig_init(); ****/
1844 	havedig = 0;
1845 	s0 = *(CONST unsigned char **)sp + 2;
1846 	while(s0[havedig] == '0')
1847 		havedig++;
1848 	s0 += havedig;
1849 	s = s0;
1850 	decpt = 0;
1851 	zret = 0;
1852 	e = 0;
1853 	if (hexdig[*s])
1854 		havedig++;
1855 	else {
1856 		zret = 1;
1857 #ifdef USE_LOCALE
1858 		for(i = 0; decimalpoint[i]; ++i) {
1859 			if (s[i] != decimalpoint[i])
1860 				goto pcheck;
1861 			}
1862 		decpt = s += i;
1863 #else
1864 		if (*s != '.')
1865 			goto pcheck;
1866 		decpt = ++s;
1867 #endif
1868 		if (!hexdig[*s])
1869 			goto pcheck;
1870 		while(*s == '0')
1871 			s++;
1872 		if (hexdig[*s])
1873 			zret = 0;
1874 		havedig = 1;
1875 		s0 = s;
1876 		}
1877 	while(hexdig[*s])
1878 		s++;
1879 #ifdef USE_LOCALE
1880 	if (*s == *decimalpoint && !decpt) {
1881 		for(i = 1; decimalpoint[i]; ++i) {
1882 			if (s[i] != decimalpoint[i])
1883 				goto pcheck;
1884 			}
1885 		decpt = s += i;
1886 #else
1887 	if (*s == '.' && !decpt) {
1888 		decpt = ++s;
1889 #endif
1890 		while(hexdig[*s])
1891 			s++;
1892 		}/*}*/
1893 	if (decpt)
1894 		e = -(((Long)(s-decpt)) << 2);
1895  pcheck:
1896 	s1 = s;
1897 	big = esign = 0;
1898 	switch(*s) {
1899 	  case 'p':
1900 	  case 'P':
1901 		switch(*++s) {
1902 		  case '-':
1903 			esign = 1;
1904 			ZEND_FALLTHROUGH;
1905 		  case '+':
1906 			s++;
1907 		  }
1908 		if ((n = hexdig[*s]) == 0 || n > 0x19) {
1909 			s = s1;
1910 			break;
1911 			}
1912 		e1 = n - 0x10;
1913 		while((n = hexdig[*++s]) !=0 && n <= 0x19) {
1914 			if (e1 & 0xf8000000)
1915 				big = 1;
1916 			e1 = 10*e1 + n - 0x10;
1917 			}
1918 		if (esign)
1919 			e1 = -e1;
1920 		e += e1;
1921 	  }
1922 	*sp = (char*)s;
1923 	if (!havedig)
1924 		*sp = (char*)s0 - 1;
1925 	if (zret)
1926 		goto retz1;
1927 	if (big) {
1928 		if (esign) {
1929 #ifdef IEEE_Arith
1930 			switch(rounding) {
1931 			  case Round_up:
1932 				if (sign)
1933 					break;
1934 				goto ret_tiny;
1935 			  case Round_down:
1936 				if (!sign)
1937 					break;
1938 				goto ret_tiny;
1939 			  }
1940 #endif
1941 			goto retz;
1942 #ifdef IEEE_Arith
1943  ret_tinyf:
1944 			Bfree(b);
1945  ret_tiny:
1946 #ifndef NO_ERRNO
1947 			errno = ERANGE;
1948 #endif
1949 			word0(rvp) = 0;
1950 			word1(rvp) = 1;
1951 			return;
1952 #endif /* IEEE_Arith */
1953 			}
1954 		switch(rounding) {
1955 		  case Round_near:
1956 			goto ovfl1;
1957 		  case Round_up:
1958 			if (!sign)
1959 				goto ovfl1;
1960 			goto ret_big;
1961 		  case Round_down:
1962 			if (sign)
1963 				goto ovfl1;
1964 			goto ret_big;
1965 		  }
1966  ret_big:
1967 		word0(rvp) = Big0;
1968 		word1(rvp) = Big1;
1969 		return;
1970 		}
1971 	n = s1 - s0 - 1;
1972 	for(k = 0; n > (1 << (kshift-2)) - 1; n >>= 1)
1973 		k++;
1974 	b = Balloc(k);
1975 	x = b->x;
1976 	n = 0;
1977 	L = 0;
1978 #ifdef USE_LOCALE
1979 	for(i = 0; decimalpoint[i+1]; ++i);
1980 #endif
1981 	while(s1 > s0) {
1982 #ifdef USE_LOCALE
1983 		if (*--s1 == decimalpoint[i]) {
1984 			s1 -= i;
1985 			continue;
1986 			}
1987 #else
1988 		if (*--s1 == '.')
1989 			continue;
1990 #endif
1991 		if (n == ULbits) {
1992 			*x++ = L;
1993 			L = 0;
1994 			n = 0;
1995 			}
1996 		L |= (hexdig[*s1] & 0x0f) << n;
1997 		n += 4;
1998 		}
1999 	*x++ = L;
2000 	b->wds = n = x - b->x;
2001 	n = ULbits*n - hi0bits(L);
2002 	nbits = Nbits;
2003 	lostbits = 0;
2004 	x = b->x;
2005 	if (n > nbits) {
2006 		n -= nbits;
2007 		if (any_on(b,n)) {
2008 			lostbits = 1;
2009 			k = n - 1;
2010 			if (x[k>>kshift] & 1 << (k & kmask)) {
2011 				lostbits = 2;
2012 				if (k > 0 && any_on(b,k))
2013 					lostbits = 3;
2014 				}
2015 			}
2016 		rshift(b, n);
2017 		e += n;
2018 		}
2019 	else if (n < nbits) {
2020 		n = nbits - n;
2021 		b = lshift(b, n);
2022 		e -= n;
2023 		x = b->x;
2024 		}
2025 	if (e > Emax) {
2026  ovfl:
2027 		Bfree(b);
2028  ovfl1:
2029 #ifndef NO_ERRNO
2030 		errno = ERANGE;
2031 #endif
2032 		word0(rvp) = Exp_mask;
2033 		word1(rvp) = 0;
2034 		return;
2035 		}
2036 	denorm = 0;
2037 	if (e < emin) {
2038 		denorm = 1;
2039 		n = emin - e;
2040 		if (n >= nbits) {
2041 #ifdef IEEE_Arith /*{*/
2042 			switch (rounding) {
2043 			  case Round_near:
2044 				if (n == nbits && (n < 2 || any_on(b,n-1)))
2045 					goto ret_tinyf;
2046 				break;
2047 			  case Round_up:
2048 				if (!sign)
2049 					goto ret_tinyf;
2050 				break;
2051 			  case Round_down:
2052 				if (sign)
2053 					goto ret_tinyf;
2054 			  }
2055 #endif /* } IEEE_Arith */
2056 			Bfree(b);
2057  retz:
2058 #ifndef NO_ERRNO
2059 			errno = ERANGE;
2060 #endif
2061  retz1:
2062 			rvp->d = 0.;
2063 			return;
2064 			}
2065 		k = n - 1;
2066 		if (lostbits)
2067 			lostbits = 1;
2068 		else if (k > 0)
2069 			lostbits = any_on(b,k);
2070 		if (x[k>>kshift] & 1 << (k & kmask))
2071 			lostbits |= 2;
2072 		nbits -= n;
2073 		rshift(b,n);
2074 		e = emin;
2075 		}
2076 	if (lostbits) {
2077 		up = 0;
2078 		switch(rounding) {
2079 		  case Round_zero:
2080 			break;
2081 		  case Round_near:
2082 			if (lostbits & 2
2083 			 && (lostbits & 1) | (x[0] & 1))
2084 				up = 1;
2085 			break;
2086 		  case Round_up:
2087 			up = 1 - sign;
2088 			break;
2089 		  case Round_down:
2090 			up = sign;
2091 		  }
2092 		if (up) {
2093 			k = b->wds;
2094 			b = increment(b);
2095 			x = b->x;
2096 			if (denorm) {
2097 #if 0
2098 				if (nbits == Nbits - 1
2099 				 && x[nbits >> kshift] & 1 << (nbits & kmask))
2100 					denorm = 0; /* not currently used */
2101 #endif
2102 				}
2103 			else if (b->wds > k
2104 			 || ((n = nbits & kmask) !=0
2105 			     && hi0bits(x[k-1]) < 32-n)) {
2106 				rshift(b,1);
2107 				if (++e > Emax)
2108 					goto ovfl;
2109 				}
2110 			}
2111 		}
2112 #ifdef IEEE_Arith
2113 	if (denorm)
2114 		word0(rvp) = b->wds > 1 ? b->x[1] & ~0x100000 : 0;
2115 	else
2116 		word0(rvp) = (b->x[1] & ~0x100000) | ((e + 0x3ff + 52) << 20);
2117 	word1(rvp) = b->x[0];
2118 #endif
2119 #ifdef IBM
2120 	if ((j = e & 3)) {
2121 		k = b->x[0] & ((1 << j) - 1);
2122 		rshift(b,j);
2123 		if (k) {
2124 			switch(rounding) {
2125 			  case Round_up:
2126 				if (!sign)
2127 					increment(b);
2128 				break;
2129 			  case Round_down:
2130 				if (sign)
2131 					increment(b);
2132 				break;
2133 			  case Round_near:
2134 				j = 1 << (j-1);
2135 				if (k & j && ((k & (j-1)) | lostbits))
2136 					increment(b);
2137 			  }
2138 			}
2139 		}
2140 	e >>= 2;
2141 	word0(rvp) = b->x[1] | ((e + 65 + 13) << 24);
2142 	word1(rvp) = b->x[0];
2143 #endif
2144 #ifdef VAX
2145 	/* The next two lines ignore swap of low- and high-order 2 bytes. */
2146 	/* word0(rvp) = (b->x[1] & ~0x800000) | ((e + 129 + 55) << 23); */
2147 	/* word1(rvp) = b->x[0]; */
2148 	word0(rvp) = ((b->x[1] & ~0x800000) >> 16) | ((e + 129 + 55) << 7) | (b->x[1] << 16);
2149 	word1(rvp) = (b->x[0] >> 16) | (b->x[0] << 16);
2150 #endif
2151 	Bfree(b);
2152 	}
2153 #endif /*!NO_HEX_FP}*/
2154 
2155  static int
2156 #ifdef KR_headers
dshift(b,p2)2157 dshift(b, p2) Bigint *b; int p2;
2158 #else
2159 dshift(Bigint *b, int p2)
2160 #endif
2161 {
2162 	int rv = hi0bits(b->x[b->wds-1]) - 4;
2163 	if (p2 > 0)
2164 		rv -= p2;
2165 	return rv & kmask;
2166 	}
2167 
2168  static int
quorem(b,S)2169 quorem
2170 #ifdef KR_headers
2171 	(b, S) Bigint *b, *S;
2172 #else
2173 	(Bigint *b, Bigint *S)
2174 #endif
2175 {
2176 	int n;
2177 	ULong *bx, *bxe, q, *sx, *sxe;
2178 #ifdef ULLong
2179 	ULLong borrow, carry, y, ys;
2180 #else
2181 	ULong borrow, carry, y, ys;
2182 #ifdef Pack_32
2183 	ULong si, z, zs;
2184 #endif
2185 #endif
2186 
2187 	n = S->wds;
2188 #ifdef DEBUG
2189 	/*debug*/ if (b->wds > n)
2190 	/*debug*/	Bug("oversize b in quorem");
2191 #endif
2192 	if (b->wds < n)
2193 		return 0;
2194 	sx = S->x;
2195 	sxe = sx + --n;
2196 	bx = b->x;
2197 	bxe = bx + n;
2198 	q = *bxe / (*sxe + 1);	/* ensure q <= true quotient */
2199 #ifdef DEBUG
2200 #ifdef NO_STRTOD_BIGCOMP
2201 	/*debug*/ if (q > 9)
2202 #else
2203 	/* An oversized q is possible when quorem is called from bigcomp and */
2204 	/* the input is near, e.g., twice the smallest denormalized number. */
2205 	/*debug*/ if (q > 15)
2206 #endif
2207 	/*debug*/	Bug("oversized quotient in quorem");
2208 #endif
2209 	if (q) {
2210 		borrow = 0;
2211 		carry = 0;
2212 		do {
2213 #ifdef ULLong
2214 			ys = *sx++ * (ULLong)q + carry;
2215 			carry = ys >> 32;
2216 			y = *bx - (ys & FFFFFFFF) - borrow;
2217 			borrow = y >> 32 & (ULong)1;
2218 			*bx++ = y & FFFFFFFF;
2219 #else
2220 #ifdef Pack_32
2221 			si = *sx++;
2222 			ys = (si & 0xffff) * q + carry;
2223 			zs = (si >> 16) * q + (ys >> 16);
2224 			carry = zs >> 16;
2225 			y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
2226 			borrow = (y & 0x10000) >> 16;
2227 			z = (*bx >> 16) - (zs & 0xffff) - borrow;
2228 			borrow = (z & 0x10000) >> 16;
2229 			Storeinc(bx, z, y);
2230 #else
2231 			ys = *sx++ * q + carry;
2232 			carry = ys >> 16;
2233 			y = *bx - (ys & 0xffff) - borrow;
2234 			borrow = (y & 0x10000) >> 16;
2235 			*bx++ = y & 0xffff;
2236 #endif
2237 #endif
2238 			}
2239 			while(sx <= sxe);
2240 		if (!*bxe) {
2241 			bx = b->x;
2242 			while(--bxe > bx && !*bxe)
2243 				--n;
2244 			b->wds = n;
2245 			}
2246 		}
2247 	if (cmp(b, S) >= 0) {
2248 		q++;
2249 		borrow = 0;
2250 		carry = 0;
2251 		bx = b->x;
2252 		sx = S->x;
2253 		do {
2254 #ifdef ULLong
2255 			ys = *sx++ + carry;
2256 			carry = ys >> 32;
2257 			y = *bx - (ys & FFFFFFFF) - borrow;
2258 			borrow = y >> 32 & (ULong)1;
2259 			*bx++ = y & FFFFFFFF;
2260 #else
2261 #ifdef Pack_32
2262 			si = *sx++;
2263 			ys = (si & 0xffff) + carry;
2264 			zs = (si >> 16) + (ys >> 16);
2265 			carry = zs >> 16;
2266 			y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
2267 			borrow = (y & 0x10000) >> 16;
2268 			z = (*bx >> 16) - (zs & 0xffff) - borrow;
2269 			borrow = (z & 0x10000) >> 16;
2270 			Storeinc(bx, z, y);
2271 #else
2272 			ys = *sx++ + carry;
2273 			carry = ys >> 16;
2274 			y = *bx - (ys & 0xffff) - borrow;
2275 			borrow = (y & 0x10000) >> 16;
2276 			*bx++ = y & 0xffff;
2277 #endif
2278 #endif
2279 			}
2280 			while(sx <= sxe);
2281 		bx = b->x;
2282 		bxe = bx + n;
2283 		if (!*bxe) {
2284 			while(--bxe > bx && !*bxe)
2285 				--n;
2286 			b->wds = n;
2287 			}
2288 		}
2289 	return q;
2290 	}
2291 
2292 #if defined(Avoid_Underflow) || !defined(NO_STRTOD_BIGCOMP) /*{*/
2293  static double
sulp(x,bc)2294 sulp
2295 #ifdef KR_headers
2296 	(x, bc) U *x; BCinfo *bc;
2297 #else
2298 	(U *x, BCinfo *bc)
2299 #endif
2300 {
2301 	U u;
2302 	double rv;
2303 	int i;
2304 
2305 	rv = ulp(x);
2306 	if (!bc->scale || (i = 2*P + 1 - ((word0(x) & Exp_mask) >> Exp_shift)) <= 0)
2307 		return rv; /* Is there an example where i <= 0 ? */
2308 	word0(&u) = Exp_1 + (i << Exp_shift);
2309 	word1(&u) = 0;
2310 	return rv * u.d;
2311 	}
2312 #endif /*}*/
2313 
2314 #ifndef NO_STRTOD_BIGCOMP
2315  static void
bigcomp(rv,s0,bc)2316 bigcomp
2317 #ifdef KR_headers
2318 	(rv, s0, bc)
2319 	U *rv; CONST char *s0; BCinfo *bc;
2320 #else
2321 	(U *rv, const char *s0, BCinfo *bc)
2322 #endif
2323 {
2324 	Bigint *b, *d;
2325 	int b2, bbits, d2, dd, dig, dsign, i, j, nd, nd0, p2, p5, speccase;
2326 
2327 	dsign = bc->dsign;
2328 	nd = bc->nd;
2329 	nd0 = bc->nd0;
2330 	p5 = nd + bc->e0 - 1;
2331 	speccase = 0;
2332 #ifndef Sudden_Underflow
2333 	if (rv->d == 0.) {	/* special case: value near underflow-to-zero */
2334 				/* threshold was rounded to zero */
2335 		b = i2b(1);
2336 		p2 = Emin - P + 1;
2337 		bbits = 1;
2338 #ifdef Avoid_Underflow
2339 		word0(rv) = (P+2) << Exp_shift;
2340 #else
2341 		word1(rv) = 1;
2342 #endif
2343 		i = 0;
2344 #ifdef Honor_FLT_ROUNDS
2345 		if (bc->rounding == 1)
2346 #endif
2347 			{
2348 			speccase = 1;
2349 			--p2;
2350 			dsign = 0;
2351 			goto have_i;
2352 			}
2353 		}
2354 	else
2355 #endif
2356 		b = d2b(rv, &p2, &bbits);
2357 #ifdef Avoid_Underflow
2358 	p2 -= bc->scale;
2359 #endif
2360 	/* floor(log2(rv)) == bbits - 1 + p2 */
2361 	/* Check for denormal case. */
2362 	i = P - bbits;
2363 	if (i > (j = P - Emin - 1 + p2)) {
2364 #ifdef Sudden_Underflow
2365 		Bfree(b);
2366 		b = i2b(1);
2367 		p2 = Emin;
2368 		i = P - 1;
2369 #ifdef Avoid_Underflow
2370 		word0(rv) = (1 + bc->scale) << Exp_shift;
2371 #else
2372 		word0(rv) = Exp_msk1;
2373 #endif
2374 		word1(rv) = 0;
2375 #else
2376 		i = j;
2377 #endif
2378 		}
2379 #ifdef Honor_FLT_ROUNDS
2380 	if (bc->rounding != 1) {
2381 		if (i > 0)
2382 			b = lshift(b, i);
2383 		if (dsign)
2384 			b = increment(b);
2385 		}
2386 	else
2387 #endif
2388 		{
2389 		b = lshift(b, ++i);
2390 		b->x[0] |= 1;
2391 		}
2392 #ifndef Sudden_Underflow
2393  have_i:
2394 #endif
2395 	p2 -= p5 + i;
2396 	d = i2b(1);
2397 	/* Arrange for convenient computation of quotients:
2398 	 * shift left if necessary so divisor has 4 leading 0 bits.
2399 	 */
2400 	if (p5 > 0)
2401 		d = pow5mult(d, p5);
2402 	else if (p5 < 0)
2403 		b = pow5mult(b, -p5);
2404 	if (p2 > 0) {
2405 		b2 = p2;
2406 		d2 = 0;
2407 		}
2408 	else {
2409 		b2 = 0;
2410 		d2 = -p2;
2411 		}
2412 	i = dshift(d, d2);
2413 	if ((b2 += i) > 0)
2414 		b = lshift(b, b2);
2415 	if ((d2 += i) > 0)
2416 		d = lshift(d, d2);
2417 
2418 	/* Now b/d = exactly half-way between the two floating-point values */
2419 	/* on either side of the input string.  Compute first digit of b/d. */
2420 
2421 	if (!(dig = quorem(b,d))) {
2422 		b = multadd(b, 10, 0);	/* very unlikely */
2423 		dig = quorem(b,d);
2424 		}
2425 
2426 	/* Compare b/d with s0 */
2427 
2428 	for(i = 0; i < nd0; ) {
2429 		if ((dd = s0[i++] - '0' - dig))
2430 			goto ret;
2431 		if (!b->x[0] && b->wds == 1) {
2432 			if (i < nd)
2433 				dd = 1;
2434 			goto ret;
2435 			}
2436 		b = multadd(b, 10, 0);
2437 		dig = quorem(b,d);
2438 		}
2439 	for(j = bc->dp1; i++ < nd;) {
2440 		if ((dd = s0[j++] - '0' - dig))
2441 			goto ret;
2442 		if (!b->x[0] && b->wds == 1) {
2443 			if (i < nd)
2444 				dd = 1;
2445 			goto ret;
2446 			}
2447 		b = multadd(b, 10, 0);
2448 		dig = quorem(b,d);
2449 		}
2450 	if (dig > 0 || b->x[0] || b->wds > 1)
2451 		dd = -1;
2452  ret:
2453 	Bfree(b);
2454 	Bfree(d);
2455 #ifdef Honor_FLT_ROUNDS
2456 	if (bc->rounding != 1) {
2457 		if (dd < 0) {
2458 			if (bc->rounding == 0) {
2459 				if (!dsign)
2460 					goto retlow1;
2461 				}
2462 			else if (dsign)
2463 				goto rethi1;
2464 			}
2465 		else if (dd > 0) {
2466 			if (bc->rounding == 0) {
2467 				if (dsign)
2468 					goto rethi1;
2469 				goto ret1;
2470 				}
2471 			if (!dsign)
2472 				goto rethi1;
2473 			dval(rv) += 2.*sulp(rv,bc);
2474 			}
2475 		else {
2476 			bc->inexact = 0;
2477 			if (dsign)
2478 				goto rethi1;
2479 			}
2480 		}
2481 	else
2482 #endif
2483 	if (speccase) {
2484 		if (dd <= 0)
2485 			rv->d = 0.;
2486 		}
2487 	else if (dd < 0) {
2488 		if (!dsign)	/* does not happen for round-near */
2489 retlow1:
2490 			dval(rv) -= sulp(rv,bc);
2491 		}
2492 	else if (dd > 0) {
2493 		if (dsign) {
2494  rethi1:
2495 			dval(rv) += sulp(rv,bc);
2496 			}
2497 		}
2498 	else {
2499 		/* Exact half-way case:  apply round-even rule. */
2500 		if ((j = ((word0(rv) & Exp_mask) >> Exp_shift) - bc->scale) <= 0) {
2501 			i = 1 - j;
2502 			if (i <= 31) {
2503 				if (word1(rv) & (0x1 << i))
2504 					goto odd;
2505 				}
2506 			else if (word0(rv) & (0x1 << (i-32)))
2507 				goto odd;
2508 			}
2509 		else if (word1(rv) & 1) {
2510  odd:
2511 			if (dsign)
2512 				goto rethi1;
2513 			goto retlow1;
2514 			}
2515 		}
2516 
2517 #ifdef Honor_FLT_ROUNDS
2518  ret1:
2519 #endif
2520 	return;
2521 	}
2522 #endif /* NO_STRTOD_BIGCOMP */
2523 
2524 ZEND_API double
zend_strtod(s00,se)2525 zend_strtod
2526 #ifdef KR_headers
2527 	(s00, se) CONST char *s00; char **se;
2528 #else
2529 	(const char *s00, const char **se)
2530 #endif
2531 {
2532 	int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, e, e1;
2533 	int esign, i, j, k, nd, nd0, nf, nz, nz0, nz1, sign;
2534 	CONST char *s, *s0, *s1;
2535 	volatile double aadj, aadj1;
2536 	Long L;
2537 	U aadj2, adj, rv, rv0;
2538 	ULong y, z;
2539 	BCinfo bc;
2540 	Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
2541 #ifdef Avoid_Underflow
2542 	ULong Lsb, Lsb1;
2543 #endif
2544 #ifdef SET_INEXACT
2545 	int oldinexact;
2546 #endif
2547 #ifndef NO_STRTOD_BIGCOMP
2548 	int req_bigcomp = 0;
2549 #endif
2550 #ifdef Honor_FLT_ROUNDS /*{*/
2551 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
2552 	bc.rounding = Flt_Rounds;
2553 #else /*}{*/
2554 	bc.rounding = 1;
2555 	switch(fegetround()) {
2556 	  case FE_TOWARDZERO:	bc.rounding = 0; break;
2557 	  case FE_UPWARD:	bc.rounding = 2; break;
2558 	  case FE_DOWNWARD:	bc.rounding = 3;
2559 	  }
2560 #endif /*}}*/
2561 #endif /*}*/
2562 #ifdef USE_LOCALE
2563 	CONST char *s2;
2564 #endif
2565 
2566 	sign = nz0 = nz1 = nz = bc.dplen = bc.uflchk = 0;
2567 	dval(&rv) = 0.;
2568 	for(s = s00;;s++) switch(*s) {
2569 		case '-':
2570 			sign = 1;
2571 			ZEND_FALLTHROUGH;
2572 		case '+':
2573 			if (*++s)
2574 				goto break2;
2575 			ZEND_FALLTHROUGH;
2576 		case 0:
2577 			goto ret0;
2578 		case '\t':
2579 		case '\n':
2580 		case '\v':
2581 		case '\f':
2582 		case '\r':
2583 		case ' ':
2584 			continue;
2585 		default:
2586 			goto break2;
2587 		}
2588  break2:
2589 	if (*s == '0') {
2590 #ifndef NO_HEX_FP /*{*/
2591 		switch(s[1]) {
2592 		  case 'x':
2593 		  case 'X':
2594 #ifdef Honor_FLT_ROUNDS
2595 			gethex(&s, &rv, bc.rounding, sign);
2596 #else
2597 			gethex(&s, &rv, 1, sign);
2598 #endif
2599 			goto ret;
2600 		  }
2601 #endif /*}*/
2602 		nz0 = 1;
2603 		while(*++s == '0') ;
2604 		if (!*s)
2605 			goto ret;
2606 		}
2607 	s0 = s;
2608 	y = z = 0;
2609 	for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
2610 		if (nd < 9)
2611 			y = 10*y + c - '0';
2612 		else if (nd < DBL_DIG + 2)
2613 			z = 10*z + c - '0';
2614 	nd0 = nd;
2615 	bc.dp0 = bc.dp1 = s - s0;
2616 	for(s1 = s; s1 > s0 && *--s1 == '0'; )
2617 		++nz1;
2618 #ifdef USE_LOCALE
2619 	s1 = localeconv()->decimal_point;
2620 	if (c == *s1) {
2621 		c = '.';
2622 		if (*++s1) {
2623 			s2 = s;
2624 			for(;;) {
2625 				if (*++s2 != *s1) {
2626 					c = 0;
2627 					break;
2628 					}
2629 				if (!*++s1) {
2630 					s = s2;
2631 					break;
2632 					}
2633 				}
2634 			}
2635 		}
2636 #endif
2637 	if (c == '.') {
2638 		c = *++s;
2639 		bc.dp1 = s - s0;
2640 		bc.dplen = bc.dp1 - bc.dp0;
2641 		if (!nd) {
2642 			for(; c == '0'; c = *++s)
2643 				nz++;
2644 			if (c > '0' && c <= '9') {
2645 				bc.dp0 = s0 - s;
2646 				bc.dp1 = bc.dp0 + bc.dplen;
2647 				s0 = s;
2648 				nf += nz;
2649 				nz = 0;
2650 				goto have_dig;
2651 				}
2652 			goto dig_done;
2653 			}
2654 		for(; c >= '0' && c <= '9'; c = *++s) {
2655  have_dig:
2656 			nz++;
2657 			if (c -= '0') {
2658 				nf += nz;
2659 				for(i = 1; i < nz; i++)
2660 					if (nd++ < 9)
2661 						y *= 10;
2662 					else if (nd <= DBL_DIG + 2)
2663 						z *= 10;
2664 				if (nd++ < 9)
2665 					y = 10*y + c;
2666 				else if (nd <= DBL_DIG + 2)
2667 					z = 10*z + c;
2668 				nz = nz1 = 0;
2669 				}
2670 			}
2671 		}
2672  dig_done:
2673  	if (nd < 0) {
2674  		/* overflow */
2675  		nd = DBL_DIG + 2;
2676  	}
2677  	if (nf < 0) {
2678  		/* overflow */
2679  		nf = DBL_DIG + 2;
2680  	}
2681 	e = 0;
2682 	if (c == 'e' || c == 'E') {
2683 		if (!nd && !nz && !nz0) {
2684 			goto ret0;
2685 			}
2686 		s00 = s;
2687 		esign = 0;
2688 		switch(c = *++s) {
2689 			case '-':
2690 				esign = 1;
2691 				ZEND_FALLTHROUGH;
2692 			case '+':
2693 				c = *++s;
2694 			}
2695 		if (c >= '0' && c <= '9') {
2696 			while(c == '0')
2697 				c = *++s;
2698 			if (c > '0' && c <= '9') {
2699 				L = c - '0';
2700 				s1 = s;
2701 				while((c = *++s) >= '0' && c <= '9')
2702 					L = (Long) (10*(ULong)L + (c - '0'));
2703 				if (s - s1 > 8 || L > 19999)
2704 					/* Avoid confusion from exponents
2705 					 * so large that e might overflow.
2706 					 */
2707 					e = 19999; /* safe for 16 bit ints */
2708 				else
2709 					e = (int)L;
2710 				if (esign)
2711 					e = -e;
2712 				}
2713 			else
2714 				e = 0;
2715 			}
2716 		else
2717 			s = s00;
2718 		}
2719 	if (!nd) {
2720 		if (!nz && !nz0) {
2721 #ifdef INFNAN_CHECK
2722 			/* Check for Nan and Infinity */
2723 			if (!bc.dplen)
2724 			 switch(c) {
2725 			  case 'i':
2726 			  case 'I':
2727 				if (match(&s,"nf")) {
2728 					--s;
2729 					if (!match(&s,"inity"))
2730 						++s;
2731 					word0(&rv) = 0x7ff00000;
2732 					word1(&rv) = 0;
2733 					goto ret;
2734 					}
2735 				break;
2736 			  case 'n':
2737 			  case 'N':
2738 				if (match(&s, "an")) {
2739 					word0(&rv) = NAN_WORD0;
2740 					word1(&rv) = NAN_WORD1;
2741 #ifndef No_Hex_NaN
2742 					if (*s == '(') /*)*/
2743 						hexnan(&rv, &s);
2744 #endif
2745 					goto ret;
2746 					}
2747 			  }
2748 #endif /* INFNAN_CHECK */
2749  ret0:
2750 			s = s00;
2751 			sign = 0;
2752 			}
2753 		goto ret;
2754 		}
2755 	bc.e0 = e1 = e -= nf;
2756 
2757 	/* Now we have nd0 digits, starting at s0, followed by a
2758 	 * decimal point, followed by nd-nd0 digits.  The number we're
2759 	 * after is the integer represented by those digits times
2760 	 * 10**e */
2761 
2762 	if (!nd0)
2763 		nd0 = nd;
2764 	k = nd < DBL_DIG + 2 ? nd : DBL_DIG + 2;
2765 	dval(&rv) = y;
2766 	if (k > 9) {
2767 #ifdef SET_INEXACT
2768 		if (k > DBL_DIG)
2769 			oldinexact = get_inexact();
2770 #endif
2771 		dval(&rv) = tens[k - 9] * dval(&rv) + z;
2772 		}
2773 	bd0 = 0;
2774 	if (nd <= DBL_DIG
2775 #ifndef RND_PRODQUOT
2776 #ifndef Honor_FLT_ROUNDS
2777 		&& Flt_Rounds == 1
2778 #endif
2779 #endif
2780 			) {
2781 		if (!e)
2782 			goto ret;
2783 #ifndef ROUND_BIASED_without_Round_Up
2784 		if (e > 0) {
2785 			if (e <= Ten_pmax) {
2786 #ifdef VAX
2787 				goto vax_ovfl_check;
2788 #else
2789 #ifdef Honor_FLT_ROUNDS
2790 				/* round correctly FLT_ROUNDS = 2 or 3 */
2791 				if (sign) {
2792 					rv.d = -rv.d;
2793 					sign = 0;
2794 					}
2795 #endif
2796 				/* rv = */ rounded_product(dval(&rv), tens[e]);
2797 				goto ret;
2798 #endif
2799 				}
2800 			i = DBL_DIG - nd;
2801 			if (e <= Ten_pmax + i) {
2802 				/* A fancier test would sometimes let us do
2803 				 * this for larger i values.
2804 				 */
2805 #ifdef Honor_FLT_ROUNDS
2806 				/* round correctly FLT_ROUNDS = 2 or 3 */
2807 				if (sign) {
2808 					rv.d = -rv.d;
2809 					sign = 0;
2810 					}
2811 #endif
2812 				e -= i;
2813 				dval(&rv) *= tens[i];
2814 #ifdef VAX
2815 				/* VAX exponent range is so narrow we must
2816 				 * worry about overflow here...
2817 				 */
2818  vax_ovfl_check:
2819 				word0(&rv) -= P*Exp_msk1;
2820 				/* rv = */ rounded_product(dval(&rv), tens[e]);
2821 				if ((word0(&rv) & Exp_mask)
2822 				 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
2823 					goto ovfl;
2824 				word0(&rv) += P*Exp_msk1;
2825 #else
2826 				/* rv = */ rounded_product(dval(&rv), tens[e]);
2827 #endif
2828 				goto ret;
2829 				}
2830 			}
2831 #ifndef Inaccurate_Divide
2832 		else if (e >= -Ten_pmax) {
2833 #ifdef Honor_FLT_ROUNDS
2834 			/* round correctly FLT_ROUNDS = 2 or 3 */
2835 			if (sign) {
2836 				rv.d = -rv.d;
2837 				sign = 0;
2838 				}
2839 #endif
2840 			/* rv = */ rounded_quotient(dval(&rv), tens[-e]);
2841 			goto ret;
2842 			}
2843 #endif
2844 #endif /* ROUND_BIASED_without_Round_Up */
2845 		}
2846 	e1 += nd - k;
2847 
2848 #ifdef IEEE_Arith
2849 #ifdef SET_INEXACT
2850 	bc.inexact = 1;
2851 	if (k <= DBL_DIG)
2852 		oldinexact = get_inexact();
2853 #endif
2854 #ifdef Avoid_Underflow
2855 	bc.scale = 0;
2856 #endif
2857 #ifdef Honor_FLT_ROUNDS
2858 	if (bc.rounding >= 2) {
2859 		if (sign)
2860 			bc.rounding = bc.rounding == 2 ? 0 : 2;
2861 		else
2862 			if (bc.rounding != 2)
2863 				bc.rounding = 0;
2864 		}
2865 #endif
2866 #endif /*IEEE_Arith*/
2867 
2868 	/* Get starting approximation = rv * 10**e1 */
2869 
2870 	if (e1 > 0) {
2871 		if ((i = e1 & 15))
2872 			dval(&rv) *= tens[i];
2873 		if (e1 &= ~15) {
2874 			if (e1 > DBL_MAX_10_EXP) {
2875  ovfl:
2876 				/* Can't trust HUGE_VAL */
2877 #ifdef IEEE_Arith
2878 #ifdef Honor_FLT_ROUNDS
2879 				switch(bc.rounding) {
2880 				  case 0: /* toward 0 */
2881 				  case 3: /* toward -infinity */
2882 					word0(&rv) = Big0;
2883 					word1(&rv) = Big1;
2884 					break;
2885 				  default:
2886 					word0(&rv) = Exp_mask;
2887 					word1(&rv) = 0;
2888 				  }
2889 #else /*Honor_FLT_ROUNDS*/
2890 				word0(&rv) = Exp_mask;
2891 				word1(&rv) = 0;
2892 #endif /*Honor_FLT_ROUNDS*/
2893 #ifdef SET_INEXACT
2894 				/* set overflow bit */
2895 				dval(&rv0) = 1e300;
2896 				dval(&rv0) *= dval(&rv0);
2897 #endif
2898 #else /*IEEE_Arith*/
2899 				word0(&rv) = Big0;
2900 				word1(&rv) = Big1;
2901 #endif /*IEEE_Arith*/
2902  range_err:
2903 				if (bd0) {
2904 					Bfree(bb);
2905 					Bfree(bd);
2906 					Bfree(bs);
2907 					Bfree(bd0);
2908 					Bfree(delta);
2909 					}
2910 #ifndef NO_ERRNO
2911 				errno = ERANGE;
2912 #endif
2913 				goto ret;
2914 				}
2915 			e1 >>= 4;
2916 			for(j = 0; e1 > 1; j++, e1 >>= 1)
2917 				if (e1 & 1)
2918 					dval(&rv) *= bigtens[j];
2919 		/* The last multiplication could overflow. */
2920 			word0(&rv) -= P*Exp_msk1;
2921 			dval(&rv) *= bigtens[j];
2922 			if ((z = word0(&rv) & Exp_mask)
2923 			 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
2924 				goto ovfl;
2925 			if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
2926 				/* set to largest number */
2927 				/* (Can't trust DBL_MAX) */
2928 				word0(&rv) = Big0;
2929 				word1(&rv) = Big1;
2930 				}
2931 			else
2932 				word0(&rv) += P*Exp_msk1;
2933 			}
2934 		}
2935 	else if (e1 < 0) {
2936 		e1 = -e1;
2937 		if ((i = e1 & 15))
2938 			dval(&rv) /= tens[i];
2939 		if (e1 >>= 4) {
2940 			if (e1 >= 1 << n_bigtens)
2941 				goto undfl;
2942 #ifdef Avoid_Underflow
2943 			if (e1 & Scale_Bit)
2944 				bc.scale = 2*P;
2945 			for(j = 0; e1 > 0; j++, e1 >>= 1)
2946 				if (e1 & 1)
2947 					dval(&rv) *= tinytens[j];
2948 			if (bc.scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask)
2949 						>> Exp_shift)) > 0) {
2950 				/* scaled rv is denormal; clear j low bits */
2951 				if (j >= 32) {
2952 					if (j > 54)
2953 						goto undfl;
2954 					word1(&rv) = 0;
2955 					if (j >= 53)
2956 					 word0(&rv) = (P+2)*Exp_msk1;
2957 					else
2958 					 word0(&rv) &= 0xffffffff << (j-32);
2959 					}
2960 				else
2961 					word1(&rv) &= 0xffffffff << j;
2962 				}
2963 #else
2964 			for(j = 0; e1 > 1; j++, e1 >>= 1)
2965 				if (e1 & 1)
2966 					dval(&rv) *= tinytens[j];
2967 			/* The last multiplication could underflow. */
2968 			dval(&rv0) = dval(&rv);
2969 			dval(&rv) *= tinytens[j];
2970 			if (!dval(&rv)) {
2971 				dval(&rv) = 2.*dval(&rv0);
2972 				dval(&rv) *= tinytens[j];
2973 #endif
2974 				if (!dval(&rv)) {
2975  undfl:
2976 					dval(&rv) = 0.;
2977 					goto range_err;
2978 					}
2979 #ifndef Avoid_Underflow
2980 				word0(&rv) = Tiny0;
2981 				word1(&rv) = Tiny1;
2982 				/* The refinement below will clean
2983 				 * this approximation up.
2984 				 */
2985 				}
2986 #endif
2987 			}
2988 		}
2989 
2990 	/* Now the hard part -- adjusting rv to the correct value.*/
2991 
2992 	/* Put digits into bd: true value = bd * 10^e */
2993 
2994 	bc.nd = nd - nz1;
2995 #ifndef NO_STRTOD_BIGCOMP
2996 	bc.nd0 = nd0;	/* Only needed if nd > strtod_diglim, but done here */
2997 			/* to silence an erroneous warning about bc.nd0 */
2998 			/* possibly not being initialized. */
2999 	if (nd > strtod_diglim) {
3000 		/* ASSERT(strtod_diglim >= 18); 18 == one more than the */
3001 		/* minimum number of decimal digits to distinguish double values */
3002 		/* in IEEE arithmetic. */
3003 		i = j = 18;
3004 		if (i > nd0)
3005 			j += bc.dplen;
3006 		for(;;) {
3007 			if (--j < bc.dp1 && j >= bc.dp0)
3008 				j = bc.dp0 - 1;
3009 			if (s0[j] != '0')
3010 				break;
3011 			--i;
3012 			}
3013 		e += nd - i;
3014 		nd = i;
3015 		if (nd0 > nd)
3016 			nd0 = nd;
3017 		if (nd < 9) { /* must recompute y */
3018 			y = 0;
3019 			for(i = 0; i < nd0; ++i)
3020 				y = 10*y + s0[i] - '0';
3021 			for(j = bc.dp1; i < nd; ++i)
3022 				y = 10*y + s0[j++] - '0';
3023 			}
3024 		}
3025 #endif
3026 	bd0 = s2b(s0, nd0, nd, y, bc.dplen);
3027 
3028 	for(;;) {
3029 		bd = Balloc(bd0->k);
3030 		Bcopy(bd, bd0);
3031 		bb = d2b(&rv, &bbe, &bbbits);	/* rv = bb * 2^bbe */
3032 		bs = i2b(1);
3033 
3034 		if (e >= 0) {
3035 			bb2 = bb5 = 0;
3036 			bd2 = bd5 = e;
3037 			}
3038 		else {
3039 			bb2 = bb5 = -e;
3040 			bd2 = bd5 = 0;
3041 			}
3042 		if (bbe >= 0)
3043 			bb2 += bbe;
3044 		else
3045 			bd2 -= bbe;
3046 		bs2 = bb2;
3047 #ifdef Honor_FLT_ROUNDS
3048 		if (bc.rounding != 1)
3049 			bs2++;
3050 #endif
3051 #ifdef Avoid_Underflow
3052 		Lsb = LSB;
3053 		Lsb1 = 0;
3054 		j = bbe - bc.scale;
3055 		i = j + bbbits - 1;	/* logb(rv) */
3056 		j = P + 1 - bbbits;
3057 		if (i < Emin) {	/* denormal */
3058 			i = Emin - i;
3059 			j -= i;
3060 			if (i < 32)
3061 				Lsb <<= i;
3062 			else if (i < 52)
3063 				Lsb1 = Lsb << (i-32);
3064 			else
3065 				Lsb1 = Exp_mask;
3066 			}
3067 #else /*Avoid_Underflow*/
3068 #ifdef Sudden_Underflow
3069 #ifdef IBM
3070 		j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
3071 #else
3072 		j = P + 1 - bbbits;
3073 #endif
3074 #else /*Sudden_Underflow*/
3075 		j = bbe;
3076 		i = j + bbbits - 1;	/* logb(rv) */
3077 		if (i < Emin)	/* denormal */
3078 			j += P - Emin;
3079 		else
3080 			j = P + 1 - bbbits;
3081 #endif /*Sudden_Underflow*/
3082 #endif /*Avoid_Underflow*/
3083 		bb2 += j;
3084 		bd2 += j;
3085 #ifdef Avoid_Underflow
3086 		bd2 += bc.scale;
3087 #endif
3088 		i = bb2 < bd2 ? bb2 : bd2;
3089 		if (i > bs2)
3090 			i = bs2;
3091 		if (i > 0) {
3092 			bb2 -= i;
3093 			bd2 -= i;
3094 			bs2 -= i;
3095 			}
3096 		if (bb5 > 0) {
3097 			bs = pow5mult(bs, bb5);
3098 			bb1 = mult(bs, bb);
3099 			Bfree(bb);
3100 			bb = bb1;
3101 			}
3102 		if (bb2 > 0)
3103 			bb = lshift(bb, bb2);
3104 		if (bd5 > 0)
3105 			bd = pow5mult(bd, bd5);
3106 		if (bd2 > 0)
3107 			bd = lshift(bd, bd2);
3108 		if (bs2 > 0)
3109 			bs = lshift(bs, bs2);
3110 		delta = diff(bb, bd);
3111 		bc.dsign = delta->sign;
3112 		delta->sign = 0;
3113 		i = cmp(delta, bs);
3114 #ifndef NO_STRTOD_BIGCOMP /*{*/
3115 		if (bc.nd > nd && i <= 0) {
3116 			if (bc.dsign) {
3117 				/* Must use bigcomp(). */
3118 				req_bigcomp = 1;
3119 				break;
3120 				}
3121 #ifdef Honor_FLT_ROUNDS
3122 			if (bc.rounding != 1) {
3123 				if (i < 0) {
3124 					req_bigcomp = 1;
3125 					break;
3126 					}
3127 				}
3128 			else
3129 #endif
3130 				i = -1;	/* Discarded digits make delta smaller. */
3131 			}
3132 #endif /*}*/
3133 #ifdef Honor_FLT_ROUNDS /*{*/
3134 		if (bc.rounding != 1) {
3135 			if (i < 0) {
3136 				/* Error is less than an ulp */
3137 				if (!delta->x[0] && delta->wds <= 1) {
3138 					/* exact */
3139 #ifdef SET_INEXACT
3140 					bc.inexact = 0;
3141 #endif
3142 					break;
3143 					}
3144 				if (bc.rounding) {
3145 					if (bc.dsign) {
3146 						adj.d = 1.;
3147 						goto apply_adj;
3148 						}
3149 					}
3150 				else if (!bc.dsign) {
3151 					adj.d = -1.;
3152 					if (!word1(&rv)
3153 					 && !(word0(&rv) & Frac_mask)) {
3154 						y = word0(&rv) & Exp_mask;
3155 #ifdef Avoid_Underflow
3156 						if (!bc.scale || y > 2*P*Exp_msk1)
3157 #else
3158 						if (y)
3159 #endif
3160 						  {
3161 						  delta = lshift(delta,Log2P);
3162 						  if (cmp(delta, bs) <= 0)
3163 							adj.d = -0.5;
3164 						  }
3165 						}
3166  apply_adj:
3167 #ifdef Avoid_Underflow /*{*/
3168 					if (bc.scale && (y = word0(&rv) & Exp_mask)
3169 						<= 2*P*Exp_msk1)
3170 					  word0(&adj) += (2*P+1)*Exp_msk1 - y;
3171 #else
3172 #ifdef Sudden_Underflow
3173 					if ((word0(&rv) & Exp_mask) <=
3174 							P*Exp_msk1) {
3175 						word0(&rv) += P*Exp_msk1;
3176 						dval(&rv) += adj.d*ulp(dval(&rv));
3177 						word0(&rv) -= P*Exp_msk1;
3178 						}
3179 					else
3180 #endif /*Sudden_Underflow*/
3181 #endif /*Avoid_Underflow}*/
3182 					dval(&rv) += adj.d*ulp(&rv);
3183 					}
3184 				break;
3185 				}
3186 			adj.d = ratio(delta, bs);
3187 			if (adj.d < 1.)
3188 				adj.d = 1.;
3189 			if (adj.d <= 0x7ffffffe) {
3190 				/* adj = rounding ? ceil(adj) : floor(adj); */
3191 				y = adj.d;
3192 				if (y != adj.d) {
3193 					if (!((bc.rounding>>1) ^ bc.dsign))
3194 						y++;
3195 					adj.d = y;
3196 					}
3197 				}
3198 #ifdef Avoid_Underflow /*{*/
3199 			if (bc.scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
3200 				word0(&adj) += (2*P+1)*Exp_msk1 - y;
3201 #else
3202 #ifdef Sudden_Underflow
3203 			if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
3204 				word0(&rv) += P*Exp_msk1;
3205 				adj.d *= ulp(dval(&rv));
3206 				if (bc.dsign)
3207 					dval(&rv) += adj.d;
3208 				else
3209 					dval(&rv) -= adj.d;
3210 				word0(&rv) -= P*Exp_msk1;
3211 				goto cont;
3212 				}
3213 #endif /*Sudden_Underflow*/
3214 #endif /*Avoid_Underflow}*/
3215 			adj.d *= ulp(&rv);
3216 			if (bc.dsign) {
3217 				if (word0(&rv) == Big0 && word1(&rv) == Big1)
3218 					goto ovfl;
3219 				dval(&rv) += adj.d;
3220 				}
3221 			else
3222 				dval(&rv) -= adj.d;
3223 			goto cont;
3224 			}
3225 #endif /*}Honor_FLT_ROUNDS*/
3226 
3227 		if (i < 0) {
3228 			/* Error is less than half an ulp -- check for
3229 			 * special case of mantissa a power of two.
3230 			 */
3231 			if (bc.dsign || word1(&rv) || word0(&rv) & Bndry_mask
3232 #ifdef IEEE_Arith /*{*/
3233 #ifdef Avoid_Underflow
3234 			 || (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1
3235 #else
3236 			 || (word0(&rv) & Exp_mask) <= Exp_msk1
3237 #endif
3238 #endif /*}*/
3239 				) {
3240 #ifdef SET_INEXACT
3241 				if (!delta->x[0] && delta->wds <= 1)
3242 					bc.inexact = 0;
3243 #endif
3244 				break;
3245 				}
3246 			if (!delta->x[0] && delta->wds <= 1) {
3247 				/* exact result */
3248 #ifdef SET_INEXACT
3249 				bc.inexact = 0;
3250 #endif
3251 				break;
3252 				}
3253 			delta = lshift(delta,Log2P);
3254 			if (cmp(delta, bs) > 0)
3255 				goto drop_down;
3256 			break;
3257 			}
3258 		if (i == 0) {
3259 			/* exactly half-way between */
3260 			if (bc.dsign) {
3261 				if ((word0(&rv) & Bndry_mask1) == Bndry_mask1
3262 				 &&  word1(&rv) == (
3263 #ifdef Avoid_Underflow
3264 			(bc.scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
3265 		? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
3266 #endif
3267 						   0xffffffff)) {
3268 					/*boundary case -- increment exponent*/
3269 					if (word0(&rv) == Big0 && word1(&rv) == Big1)
3270 						goto ovfl;
3271 					word0(&rv) = (word0(&rv) & Exp_mask)
3272 						+ Exp_msk1
3273 #ifdef IBM
3274 						| Exp_msk1 >> 4
3275 #endif
3276 						;
3277 					word1(&rv) = 0;
3278 #ifdef Avoid_Underflow
3279 					bc.dsign = 0;
3280 #endif
3281 					break;
3282 					}
3283 				}
3284 			else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) {
3285  drop_down:
3286 				/* boundary case -- decrement exponent */
3287 #ifdef Sudden_Underflow /*{{*/
3288 				L = word0(&rv) & Exp_mask;
3289 #ifdef IBM
3290 				if (L <  Exp_msk1)
3291 #else
3292 #ifdef Avoid_Underflow
3293 				if (L <= (bc.scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
3294 #else
3295 				if (L <= Exp_msk1)
3296 #endif /*Avoid_Underflow*/
3297 #endif /*IBM*/
3298 					{
3299 					if (bc.nd >nd) {
3300 						bc.uflchk = 1;
3301 						break;
3302 						}
3303 					goto undfl;
3304 					}
3305 				L -= Exp_msk1;
3306 #else /*Sudden_Underflow}{*/
3307 #ifdef Avoid_Underflow
3308 				if (bc.scale) {
3309 					L = word0(&rv) & Exp_mask;
3310 					if (L <= (2*P+1)*Exp_msk1) {
3311 						if (L > (P+2)*Exp_msk1)
3312 							/* round even ==> */
3313 							/* accept rv */
3314 							break;
3315 						/* rv = smallest denormal */
3316 						if (bc.nd >nd) {
3317 							bc.uflchk = 1;
3318 							break;
3319 							}
3320 						goto undfl;
3321 						}
3322 					}
3323 #endif /*Avoid_Underflow*/
3324 				L = (word0(&rv) & Exp_mask) - Exp_msk1;
3325 #endif /*Sudden_Underflow}}*/
3326 				word0(&rv) = L | Bndry_mask1;
3327 				word1(&rv) = 0xffffffff;
3328 #ifdef IBM
3329 				goto cont;
3330 #else
3331 #ifndef NO_STRTOD_BIGCOMP
3332 				if (bc.nd > nd)
3333 					goto cont;
3334 #endif
3335 				break;
3336 #endif
3337 				}
3338 #ifndef ROUND_BIASED
3339 #ifdef Avoid_Underflow
3340 			if (Lsb1) {
3341 				if (!(word0(&rv) & Lsb1))
3342 					break;
3343 				}
3344 			else if (!(word1(&rv) & Lsb))
3345 				break;
3346 #else
3347 			if (!(word1(&rv) & LSB))
3348 				break;
3349 #endif
3350 #endif
3351 			if (bc.dsign)
3352 #ifdef Avoid_Underflow
3353 				dval(&rv) += sulp(&rv, &bc);
3354 #else
3355 				dval(&rv) += ulp(&rv);
3356 #endif
3357 #ifndef ROUND_BIASED
3358 			else {
3359 #ifdef Avoid_Underflow
3360 				dval(&rv) -= sulp(&rv, &bc);
3361 #else
3362 				dval(&rv) -= ulp(&rv);
3363 #endif
3364 #ifndef Sudden_Underflow
3365 				if (!dval(&rv)) {
3366 					if (bc.nd >nd) {
3367 						bc.uflchk = 1;
3368 						break;
3369 						}
3370 					goto undfl;
3371 					}
3372 #endif
3373 				}
3374 #ifdef Avoid_Underflow
3375 			bc.dsign = 1 - bc.dsign;
3376 #endif
3377 #endif
3378 			break;
3379 			}
3380 		if ((aadj = ratio(delta, bs)) <= 2.) {
3381 			if (bc.dsign)
3382 				aadj = aadj1 = 1.;
3383 			else if (word1(&rv) || word0(&rv) & Bndry_mask) {
3384 #ifndef Sudden_Underflow
3385 				if (word1(&rv) == Tiny1 && !word0(&rv)) {
3386 					if (bc.nd >nd) {
3387 						bc.uflchk = 1;
3388 						break;
3389 						}
3390 					goto undfl;
3391 					}
3392 #endif
3393 				aadj = 1.;
3394 				aadj1 = -1.;
3395 				}
3396 			else {
3397 				/* special case -- power of FLT_RADIX to be */
3398 				/* rounded down... */
3399 
3400 				if (aadj < 2./FLT_RADIX)
3401 					aadj = 1./FLT_RADIX;
3402 				else
3403 					aadj *= 0.5;
3404 				aadj1 = -aadj;
3405 				}
3406 			}
3407 		else {
3408 			aadj *= 0.5;
3409 			aadj1 = bc.dsign ? aadj : -aadj;
3410 #ifdef Check_FLT_ROUNDS
3411 			switch(bc.rounding) {
3412 				case 2: /* towards +infinity */
3413 					aadj1 -= 0.5;
3414 					break;
3415 				case 0: /* towards 0 */
3416 				case 3: /* towards -infinity */
3417 					aadj1 += 0.5;
3418 				}
3419 #else
3420 			if (Flt_Rounds == 0)
3421 				aadj1 += 0.5;
3422 #endif /*Check_FLT_ROUNDS*/
3423 			}
3424 		y = word0(&rv) & Exp_mask;
3425 
3426 		/* Check for overflow */
3427 
3428 		if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
3429 			dval(&rv0) = dval(&rv);
3430 			word0(&rv) -= P*Exp_msk1;
3431 			adj.d = aadj1 * ulp(&rv);
3432 			dval(&rv) += adj.d;
3433 			if ((word0(&rv) & Exp_mask) >=
3434 					Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
3435 				if (word0(&rv0) == Big0 && word1(&rv0) == Big1)
3436 					goto ovfl;
3437 				word0(&rv) = Big0;
3438 				word1(&rv) = Big1;
3439 				goto cont;
3440 				}
3441 			else
3442 				word0(&rv) += P*Exp_msk1;
3443 			}
3444 		else {
3445 #ifdef Avoid_Underflow
3446 			if (bc.scale && y <= 2*P*Exp_msk1) {
3447 				if (aadj <= 0x7fffffff) {
3448 					if ((z = aadj) <= 0)
3449 						z = 1;
3450 					aadj = z;
3451 					aadj1 = bc.dsign ? aadj : -aadj;
3452 					}
3453 				dval(&aadj2) = aadj1;
3454 				word0(&aadj2) += (2*P+1)*Exp_msk1 - y;
3455 				aadj1 = dval(&aadj2);
3456 				adj.d = aadj1 * ulp(&rv);
3457 				dval(&rv) += adj.d;
3458 				if (rv.d == 0.)
3459 #ifdef NO_STRTOD_BIGCOMP
3460 					goto undfl;
3461 #else
3462 					{
3463 					req_bigcomp = 1;
3464 					break;
3465 					}
3466 #endif
3467 				}
3468 			else {
3469 				adj.d = aadj1 * ulp(&rv);
3470 				dval(&rv) += adj.d;
3471 				}
3472 #else
3473 #ifdef Sudden_Underflow
3474 			if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
3475 				dval(&rv0) = dval(&rv);
3476 				word0(&rv) += P*Exp_msk1;
3477 				adj.d = aadj1 * ulp(&rv);
3478 				dval(&rv) += adj.d;
3479 #ifdef IBM
3480 				if ((word0(&rv) & Exp_mask) <  P*Exp_msk1)
3481 #else
3482 				if ((word0(&rv) & Exp_mask) <= P*Exp_msk1)
3483 #endif
3484 					{
3485 					if (word0(&rv0) == Tiny0
3486 					 && word1(&rv0) == Tiny1) {
3487 						if (bc.nd >nd) {
3488 							bc.uflchk = 1;
3489 							break;
3490 							}
3491 						goto undfl;
3492 						}
3493 					word0(&rv) = Tiny0;
3494 					word1(&rv) = Tiny1;
3495 					goto cont;
3496 					}
3497 				else
3498 					word0(&rv) -= P*Exp_msk1;
3499 				}
3500 			else {
3501 				adj.d = aadj1 * ulp(&rv);
3502 				dval(&rv) += adj.d;
3503 				}
3504 #else /*Sudden_Underflow*/
3505 			/* Compute adj so that the IEEE rounding rules will
3506 			 * correctly round rv + adj in some half-way cases.
3507 			 * If rv * ulp(rv) is denormalized (i.e.,
3508 			 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
3509 			 * trouble from bits lost to denormalization;
3510 			 * example: 1.2e-307 .
3511 			 */
3512 			if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
3513 				aadj1 = (double)(int)(aadj + 0.5);
3514 				if (!bc.dsign)
3515 					aadj1 = -aadj1;
3516 				}
3517 			adj.d = aadj1 * ulp(&rv);
3518 			dval(&rv) += adj.d;
3519 #endif /*Sudden_Underflow*/
3520 #endif /*Avoid_Underflow*/
3521 			}
3522 		z = word0(&rv) & Exp_mask;
3523 #ifndef SET_INEXACT
3524 		if (bc.nd == nd) {
3525 #ifdef Avoid_Underflow
3526 		if (!bc.scale)
3527 #endif
3528 		if (y == z) {
3529 			/* Can we stop now? */
3530 			L = (Long)aadj;
3531 			aadj -= L;
3532 			/* The tolerances below are conservative. */
3533 			if (bc.dsign || word1(&rv) || word0(&rv) & Bndry_mask) {
3534 				if (aadj < .4999999 || aadj > .5000001)
3535 					break;
3536 				}
3537 			else if (aadj < .4999999/FLT_RADIX)
3538 				break;
3539 			}
3540 		}
3541 #endif
3542  cont:
3543 		Bfree(bb);
3544 		Bfree(bd);
3545 		Bfree(bs);
3546 		Bfree(delta);
3547 		}
3548 	Bfree(bb);
3549 	Bfree(bd);
3550 	Bfree(bs);
3551 	Bfree(bd0);
3552 	Bfree(delta);
3553 #ifndef NO_STRTOD_BIGCOMP
3554 	if (req_bigcomp) {
3555 		bd0 = 0;
3556 		bc.e0 += nz1;
3557 		bigcomp(&rv, s0, &bc);
3558 		y = word0(&rv) & Exp_mask;
3559 		if (y == Exp_mask)
3560 			goto ovfl;
3561 		if (y == 0 && rv.d == 0.)
3562 			goto undfl;
3563 		}
3564 #endif
3565 #ifdef SET_INEXACT
3566 	if (bc.inexact) {
3567 		if (!oldinexact) {
3568 			word0(&rv0) = Exp_1 + (70 << Exp_shift);
3569 			word1(&rv0) = 0;
3570 			dval(&rv0) += 1.;
3571 			}
3572 		}
3573 	else if (!oldinexact)
3574 		clear_inexact();
3575 #endif
3576 #ifdef Avoid_Underflow
3577 	if (bc.scale) {
3578 		word0(&rv0) = Exp_1 - 2*P*Exp_msk1;
3579 		word1(&rv0) = 0;
3580 		dval(&rv) *= dval(&rv0);
3581 #ifndef NO_ERRNO
3582 		/* try to avoid the bug of testing an 8087 register value */
3583 #ifdef IEEE_Arith
3584 		if (!(word0(&rv) & Exp_mask))
3585 #else
3586 		if (word0(&rv) == 0 && word1(&rv) == 0)
3587 #endif
3588 			errno = ERANGE;
3589 #endif
3590 		}
3591 #endif /* Avoid_Underflow */
3592 #ifdef SET_INEXACT
3593 	if (bc.inexact && !(word0(&rv) & Exp_mask)) {
3594 		/* set underflow bit */
3595 		dval(&rv0) = 1e-300;
3596 		dval(&rv0) *= dval(&rv0);
3597 		}
3598 #endif
3599  ret:
3600 	if (se)
3601 		*se = (char *)s;
3602 	return sign ? -dval(&rv) : dval(&rv);
3603 	}
3604 
3605 #ifndef MULTIPLE_THREADS
3606  ZEND_TLS char *dtoa_result;
3607 #endif
3608 
3609  static char *
3610 #ifdef KR_headers
rv_alloc(i)3611 rv_alloc(i) int i;
3612 #else
3613 rv_alloc(int i)
3614 #endif
3615 {
3616 	int j, k, *r;
3617 
3618 	j = sizeof(ULong);
3619 	for(k = 0;
3620 		sizeof(Bigint) - sizeof(ULong) - sizeof(int) + (size_t)j <= (size_t)i;
3621 		j <<= 1)
3622 			k++;
3623 	r = (int*)Balloc(k);
3624 	*r = k;
3625 	return
3626 #ifndef MULTIPLE_THREADS
3627 	dtoa_result =
3628 #endif
3629 		(char *)(r+1);
3630 	}
3631 
3632  static char *
3633 #ifdef KR_headers
nrv_alloc(s,rve,n)3634 nrv_alloc(s, rve, n) char *s, **rve; int n;
3635 #else
3636 nrv_alloc(const char *s, char **rve, int n)
3637 #endif
3638 {
3639 	char *rv, *t;
3640 
3641 	t = rv = rv_alloc(n);
3642 	while((*t = *s++)) t++;
3643 	if (rve)
3644 		*rve = t;
3645 	return rv;
3646 	}
3647 
3648 /* freedtoa(s) must be used to free values s returned by dtoa
3649  * when MULTIPLE_THREADS is #defined.  It should be used in all cases,
3650  * but for consistency with earlier versions of dtoa, it is optional
3651  * when MULTIPLE_THREADS is not defined.
3652  */
3653 
3654 ZEND_API void
3655 #ifdef KR_headers
zend_freedtoa(s)3656 zend_freedtoa(s) char *s;
3657 #else
3658 zend_freedtoa(char *s)
3659 #endif
3660 {
3661 	Bigint *b = (Bigint *)((int *)s - 1);
3662 	b->maxwds = 1 << (b->k = *(int*)b);
3663 	Bfree(b);
3664 #ifndef MULTIPLE_THREADS
3665 	if (s == dtoa_result)
3666 		dtoa_result = 0;
3667 #endif
3668 	}
3669 
3670 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
3671  *
3672  * Inspired by "How to Print Floating-Point Numbers Accurately" by
3673  * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
3674  *
3675  * Modifications:
3676  *	1. Rather than iterating, we use a simple numeric overestimate
3677  *	   to determine k = floor(log10(d)).  We scale relevant
3678  *	   quantities using O(log2(k)) rather than O(k) multiplications.
3679  *	2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
3680  *	   try to generate digits strictly left to right.  Instead, we
3681  *	   compute with fewer bits and propagate the carry if necessary
3682  *	   when rounding the final digit up.  This is often faster.
3683  *	3. Under the assumption that input will be rounded nearest,
3684  *	   mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
3685  *	   That is, we allow equality in stopping tests when the
3686  *	   round-nearest rule will give the same floating-point value
3687  *	   as would satisfaction of the stopping test with strict
3688  *	   inequality.
3689  *	4. We remove common factors of powers of 2 from relevant
3690  *	   quantities.
3691  *	5. When converting floating-point integers less than 1e16,
3692  *	   we use floating-point arithmetic rather than resorting
3693  *	   to multiple-precision integers.
3694  *	6. When asked to produce fewer than 15 digits, we first try
3695  *	   to get by with floating-point arithmetic; we resort to
3696  *	   multiple-precision integer arithmetic only if we cannot
3697  *	   guarantee that the floating-point calculation has given
3698  *	   the correctly rounded result.  For k requested digits and
3699  *	   "uniformly" distributed input, the probability is
3700  *	   something like 10^(k-15) that we must resort to the Long
3701  *	   calculation.
3702  */
3703 
zend_dtoa(double dd,int mode,int ndigits,int * decpt,bool * sign,char ** rve)3704 ZEND_API char *zend_dtoa(double dd, int mode, int ndigits, int *decpt, bool *sign, char **rve)
3705 {
3706  /*	Arguments ndigits, decpt, sign are similar to those
3707 	of ecvt and fcvt; trailing zeros are suppressed from
3708 	the returned string.  If not null, *rve is set to point
3709 	to the end of the return value.  If d is +-Infinity or NaN,
3710 	then *decpt is set to 9999.
3711 
3712 	mode:
3713 		0 ==> shortest string that yields d when read in
3714 			and rounded to nearest.
3715 		1 ==> like 0, but with Steele & White stopping rule;
3716 			e.g. with IEEE P754 arithmetic , mode 0 gives
3717 			1e23 whereas mode 1 gives 9.999999999999999e22.
3718 		2 ==> max(1,ndigits) significant digits.  This gives a
3719 			return value similar to that of ecvt, except
3720 			that trailing zeros are suppressed.
3721 		3 ==> through ndigits past the decimal point.  This
3722 			gives a return value similar to that from fcvt,
3723 			except that trailing zeros are suppressed, and
3724 			ndigits can be negative.
3725 		4,5 ==> similar to 2 and 3, respectively, but (in
3726 			round-nearest mode) with the tests of mode 0 to
3727 			possibly return a shorter string that rounds to d.
3728 			With IEEE arithmetic and compilation with
3729 			-DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
3730 			as modes 2 and 3 when FLT_ROUNDS != 1.
3731 		6-9 ==> Debugging modes similar to mode - 4:  don't try
3732 			fast floating-point estimate (if applicable).
3733 
3734 		Values of mode other than 0-9 are treated as mode 0.
3735 
3736 		Sufficient space is allocated to the return value
3737 		to hold the suppressed trailing zeros.
3738 	*/
3739 
3740 	int bbits, b2, b5, be, dig, i, ieps, ilim = 0, ilim0, ilim1,
3741 		j, j1 = 0, k, k0, k_check, leftright, m2, m5, s2, s5,
3742 		spec_case = 0, try_quick;
3743 	Long L;
3744 #ifndef Sudden_Underflow
3745 	int denorm;
3746 	ULong x;
3747 #endif
3748 	Bigint *b, *b1, *delta, *mlo, *mhi, *S;
3749 	U d2, eps, u;
3750 	double ds;
3751 	char *s, *s0;
3752 #ifndef No_leftright
3753 #ifdef IEEE_Arith
3754 	U eps1;
3755 #endif
3756 #endif
3757 #ifdef SET_INEXACT
3758 	int inexact, oldinexact;
3759 #endif
3760 #ifdef Honor_FLT_ROUNDS /*{*/
3761 	int Rounding;
3762 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
3763 	Rounding = Flt_Rounds;
3764 #else /*}{*/
3765 	Rounding = 1;
3766 	switch(fegetround()) {
3767 	  case FE_TOWARDZERO:	Rounding = 0; break;
3768 	  case FE_UPWARD:	Rounding = 2; break;
3769 	  case FE_DOWNWARD:	Rounding = 3;
3770 	  }
3771 #endif /*}}*/
3772 #endif /*}*/
3773 
3774 #ifndef MULTIPLE_THREADS
3775 	if (dtoa_result) {
3776