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