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