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