1 /****************************************************************
2 *
3 * The author of this software is David M. Gay.
4 *
5 * Copyright (c) 1991, 2000, 2001 by Lucent Technologies.
6 *
7 * Permission to use, copy, modify, and distribute this software for any
8 * purpose without fee is hereby granted, provided that this entire notice
9 * is included in all copies of any software which is or includes a copy
10 * or modification of this software and in all copies of the supporting
11 * documentation for such software.
12 *
13 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
14 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY
15 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
16 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
17 *
18 ***************************************************************/
19
20 /* Please send bug reports to David M. Gay (dmg at acm dot org,
21 * with " at " changed at "@" and " dot " changed to "."). */
22
23 /* On a machine with IEEE extended-precision registers, it is
24 * necessary to specify double-precision (53-bit) rounding precision
25 * before invoking strtod or dtoa. If the machine uses (the equivalent
26 * of) Intel 80x87 arithmetic, the call
27 * _control87(PC_53, MCW_PC);
28 * does this with many compilers. Whether this or another call is
29 * appropriate depends on the compiler; for this to work, it may be
30 * necessary to #include "float.h" or another system-dependent header
31 * file.
32 */
33
34 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
35 * (Note that IEEE arithmetic is disabled by gcc's -ffast-math flag.)
36 *
37 * This strtod returns a nearest machine number to the input decimal
38 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are
39 * broken by the IEEE round-even rule. Otherwise ties are broken by
40 * biased rounding (add half and chop).
41 *
42 * Inspired loosely by William D. Clinger's paper "How to Read Floating
43 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
44 *
45 * Modifications:
46 *
47 * 1. We only require IEEE, IBM, or VAX double-precision
48 * arithmetic (not IEEE double-extended).
49 * 2. We get by with floating-point arithmetic in a case that
50 * Clinger missed -- when we're computing d * 10^n
51 * for a small integer d and the integer n is not too
52 * much larger than 22 (the maximum integer k for which
53 * we can represent 10^k exactly), we may be able to
54 * compute (d*10^k) * 10^(e-k) with just one roundoff.
55 * 3. Rather than a bit-at-a-time adjustment of the binary
56 * result in the hard case, we use floating-point
57 * arithmetic to determine the adjustment to within
58 * one bit; only in really hard cases do we need to
59 * compute a second residual.
60 * 4. Because of 3., we don't need a large table of powers of 10
61 * for ten-to-e (just some small tables, e.g. of 10^k
62 * for 0 <= k <= 22).
63 */
64
65 /*
66 * #define IEEE_8087 for IEEE-arithmetic machines where the least
67 * significant byte has the lowest address.
68 * #define IEEE_MC68k for IEEE-arithmetic machines where the most
69 * significant byte has the lowest address.
70 * #define Long int on machines with 32-bit ints and 64-bit longs.
71 * #define IBM for IBM mainframe-style floating-point arithmetic.
72 * #define VAX for VAX-style floating-point arithmetic (D_floating).
73 * #define No_leftright to omit left-right logic in fast floating-point
74 * computation of dtoa. This will cause dtoa modes 4 and 5 to be
75 * treated the same as modes 2 and 3 for some inputs.
76 * #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
77 * and strtod and dtoa should round accordingly. Unless Trust_FLT_ROUNDS
78 * is also #defined, fegetround() will be queried for the rounding mode.
79 * Note that both FLT_ROUNDS and fegetround() are specified by the C99
80 * standard (and are specified to be consistent, with fesetround()
81 * affecting the value of FLT_ROUNDS), but that some (Linux) systems
82 * do not work correctly in this regard, so using fegetround() is more
83 * portable than using FLT_ROUNDS directly.
84 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
85 * and Honor_FLT_ROUNDS is not #defined.
86 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
87 * that use extended-precision instructions to compute rounded
88 * products and quotients) with IBM.
89 * #define ROUND_BIASED for IEEE-format with biased rounding and arithmetic
90 * that rounds toward +Infinity.
91 * #define ROUND_BIASED_without_Round_Up for IEEE-format with biased
92 * rounding when the underlying floating-point arithmetic uses
93 * unbiased rounding. This prevent using ordinary floating-point
94 * arithmetic when the result could be computed with one rounding error.
95 * #define Inaccurate_Divide for IEEE-format with correctly rounded
96 * products but inaccurate quotients, e.g., for Intel i860.
97 * #define NO_LONG_LONG on machines that do not have a "long long"
98 * integer type (of >= 64 bits). On such machines, you can
99 * #define Just_16 to store 16 bits per 32-bit Long when doing
100 * high-precision integer arithmetic. Whether this speeds things
101 * up or slows things down depends on the machine and the number
102 * being converted. If long long is available and the name is
103 * something other than "long long", #define Llong to be the name,
104 * and if "unsigned Llong" does not work as an unsigned version of
105 * Llong, #define #ULLong to be the corresponding unsigned type.
106 * #define KR_headers for old-style C function headers.
107 * #define Bad_float_h if your system lacks a float.h or if it does not
108 * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
109 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
110 * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
111 * if memory is available and otherwise does something you deem
112 * appropriate. If MALLOC is undefined, malloc will be invoked
113 * directly -- and assumed always to succeed. Similarly, if you
114 * want something other than the system's free() to be called to
115 * recycle memory acquired from MALLOC, #define FREE to be the
116 * name of the alternate routine. (FREE or free is only called in
117 * pathological cases, e.g., in a dtoa call after a dtoa return in
118 * mode 3 with thousands of digits requested.)
119 * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making
120 * memory allocations from a private pool of memory when possible.
121 * When used, the private pool is PRIVATE_MEM bytes long: 2304 bytes,
122 * unless #defined to be a different length. This default length
123 * suffices to get rid of MALLOC calls except for unusual cases,
124 * such as decimal-to-binary conversion of a very long string of
125 * digits. The longest string dtoa can return is about 751 bytes
126 * long. For conversions by strtod of strings of 800 digits and
127 * all dtoa conversions in single-threaded executions with 8-byte
128 * pointers, PRIVATE_MEM >= 7400 appears to suffice; with 4-byte
129 * pointers, PRIVATE_MEM >= 7112 appears adequate.
130 * #define NO_INFNAN_CHECK if you do not wish to have INFNAN_CHECK
131 * #defined automatically on IEEE systems. On such systems,
132 * when INFNAN_CHECK is #defined, strtod checks
133 * for Infinity and NaN (case insensitively). On some systems
134 * (e.g., some HP systems), it may be necessary to #define NAN_WORD0
135 * appropriately -- to the most significant word of a quiet NaN.
136 * (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.)
137 * When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined,
138 * strtod also accepts (case insensitively) strings of the form
139 * NaN(x), where x is a string of hexadecimal digits and spaces;
140 * if there is only one string of hexadecimal digits, it is taken
141 * for the 52 fraction bits of the resulting NaN; if there are two
142 * or more strings of hex digits, the first is for the high 20 bits,
143 * the second and subsequent for the low 32 bits, with intervening
144 * white space ignored; but if this results in none of the 52
145 * fraction bits being on (an IEEE Infinity symbol), then NAN_WORD0
146 * and NAN_WORD1 are used instead.
147 * #define MULTIPLE_THREADS if the system offers preemptively scheduled
148 * multiple threads. In this case, you must provide (or suitably
149 * #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed
150 * by FREE_DTOA_LOCK(n) for n = 0 or 1. (The second lock, accessed
151 * in pow5mult, ensures lazy evaluation of only one copy of high
152 * powers of 5; omitting this lock would introduce a small
153 * probability of wasting memory, but would otherwise be harmless.)
154 * You must also invoke freedtoa(s) to free the value s returned by
155 * dtoa. You may do so whether or not MULTIPLE_THREADS is #defined.
156 * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that
157 * avoids underflows on inputs whose result does not underflow.
158 * If you #define NO_IEEE_Scale on a machine that uses IEEE-format
159 * floating-point numbers and flushes underflows to zero rather
160 * than implementing gradual underflow, then you must also #define
161 * Sudden_Underflow.
162 * #define USE_LOCALE to use the current locale's decimal_point value.
163 * #define SET_INEXACT if IEEE arithmetic is being used and extra
164 * computation should be done to set the inexact flag when the
165 * result is inexact and avoid setting inexact when the result
166 * is exact. In this case, dtoa.c must be compiled in
167 * an environment, perhaps provided by #include "dtoa.c" in a
168 * suitable wrapper, that defines two functions,
169 * int get_inexact(void);
170 * void clear_inexact(void);
171 * such that get_inexact() returns a nonzero value if the
172 * inexact bit is already set, and clear_inexact() sets the
173 * inexact bit to 0. When SET_INEXACT is #defined, strtod
174 * also does extra computations to set the underflow and overflow
175 * flags when appropriate (i.e., when the result is tiny and
176 * inexact or when it is a numeric value rounded to +-infinity).
177 * #define NO_ERRNO if strtod should not assign errno = ERANGE when
178 * the result overflows to +-Infinity or underflows to 0.
179 * #define NO_HEX_FP to omit recognition of hexadecimal floating-point
180 * values by strtod.
181 * #define NO_STRTOD_BIGCOMP (on IEEE-arithmetic systems only for now)
182 * to disable logic for "fast" testing of very long input strings
183 * to strtod. This testing proceeds by initially truncating the
184 * input string, then if necessary comparing the whole string with
185 * a decimal expansion to decide close cases. This logic is only
186 * used for input more than STRTOD_DIGLIM digits long (default 40).
187 */
188
189 #include <zend_operators.h>
190 #include <zend_strtod.h>
191 #include "zend_strtod_int.h"
192
193 #ifndef Long
194 #define Long int32_t
195 #endif
196 #ifndef ULong
197 #define ULong uint32_t
198 #endif
199
200 #ifdef DEBUG
Bug(const char * message)201 static void Bug(const char *message) {
202 fprintf(stderr, "%s\n", message);
203 }
204 #endif
205
206 #include "stdlib.h"
207 #include "string.h"
208
209 #ifdef USE_LOCALE
210 #include "locale.h"
211 #endif
212
213 #ifdef Honor_FLT_ROUNDS
214 #ifndef Trust_FLT_ROUNDS
215 #include <fenv.h>
216 #endif
217 #endif
218
219 #ifdef MALLOC
220 #ifdef KR_headers
221 extern char *MALLOC();
222 #else
223 extern void *MALLOC(size_t);
224 #endif
225 #else
226 #define MALLOC malloc
227 #endif
228
229 #ifndef Omit_Private_Memory
230 #ifndef PRIVATE_MEM
231 #define PRIVATE_MEM 2304
232 #endif
233 #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
234 static double private_mem[PRIVATE_mem], *pmem_next = private_mem;
235 #endif
236
237 #undef IEEE_Arith
238 #undef Avoid_Underflow
239 #ifdef IEEE_MC68k
240 #define IEEE_Arith
241 #endif
242 #ifdef IEEE_8087
243 #define IEEE_Arith
244 #endif
245
246 #ifdef IEEE_Arith
247 #ifndef NO_INFNAN_CHECK
248 #undef INFNAN_CHECK
249 #define INFNAN_CHECK
250 #endif
251 #else
252 #undef INFNAN_CHECK
253 #define NO_STRTOD_BIGCOMP
254 #endif
255
256 #include "errno.h"
257
258 #ifdef Bad_float_h
259
260 #ifdef IEEE_Arith
261 #define DBL_DIG 15
262 #define DBL_MAX_10_EXP 308
263 #define DBL_MAX_EXP 1024
264 #define FLT_RADIX 2
265 #endif /*IEEE_Arith*/
266
267 #ifdef IBM
268 #define DBL_DIG 16
269 #define DBL_MAX_10_EXP 75
270 #define DBL_MAX_EXP 63
271 #define FLT_RADIX 16
272 #define DBL_MAX 7.2370055773322621e+75
273 #endif
274
275 #ifdef VAX
276 #define DBL_DIG 16
277 #define DBL_MAX_10_EXP 38
278 #define DBL_MAX_EXP 127
279 #define FLT_RADIX 2
280 #define DBL_MAX 1.7014118346046923e+38
281 #endif
282
283 #ifndef LONG_MAX
284 #define LONG_MAX 2147483647
285 #endif
286
287 #else /* ifndef Bad_float_h */
288 #include "float.h"
289 #endif /* Bad_float_h */
290
291 #ifndef __MATH_H__
292 #include "math.h"
293 #endif
294
295 #ifdef __cplusplus
296 extern "C" {
297 #endif
298
299 #ifndef CONST
300 #ifdef KR_headers
301 #define CONST /* blank */
302 #else
303 #define CONST const
304 #endif
305 #endif
306
307 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
308 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
309 #endif
310
311 typedef union { double d; ULong L[2]; } U;
312
313 #ifdef IEEE_8087
314 #define word0(x) (x)->L[1]
315 #define word1(x) (x)->L[0]
316 #else
317 #define word0(x) (x)->L[0]
318 #define word1(x) (x)->L[1]
319 #endif
320 #define dval(x) (x)->d
321
322 #ifndef STRTOD_DIGLIM
323 #define STRTOD_DIGLIM 40
324 #endif
325
326 #ifdef DIGLIM_DEBUG
327 extern int strtod_diglim;
328 #else
329 #define strtod_diglim STRTOD_DIGLIM
330 #endif
331
332 /* The following definition of Storeinc is appropriate for MIPS processors.
333 * An alternative that might be better on some machines is
334 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
335 */
336 #if defined(IEEE_8087) + defined(VAX) + defined(__arm__)
337 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
338 ((unsigned short *)a)[0] = (unsigned short)c, a++)
339 #else
340 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
341 ((unsigned short *)a)[1] = (unsigned short)c, a++)
342 #endif
343
344 /* #define P DBL_MANT_DIG */
345 /* Ten_pmax = floor(P*log(2)/log(5)) */
346 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
347 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
348 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
349
350 #ifdef IEEE_Arith
351 #define Exp_shift 20
352 #define Exp_shift1 20
353 #define Exp_msk1 0x100000
354 #define Exp_msk11 0x100000
355 #define Exp_mask 0x7ff00000
356 #define P 53
357 #define Nbits 53
358 #define Bias 1023
359 #define Emax 1023
360 #define Emin (-1022)
361 #define Exp_1 0x3ff00000
362 #define Exp_11 0x3ff00000
363 #define Ebits 11
364 #define Frac_mask 0xfffff
365 #define Frac_mask1 0xfffff
366 #define Ten_pmax 22
367 #define Bletch 0x10
368 #define Bndry_mask 0xfffff
369 #define Bndry_mask1 0xfffff
370 #define LSB 1
371 #define Sign_bit 0x80000000
372 #define Log2P 1
373 #define Tiny0 0
374 #define Tiny1 1
375 #define Quick_max 14
376 #define Int_max 14
377 #ifndef NO_IEEE_Scale
378 #define Avoid_Underflow
379 #ifdef Flush_Denorm /* debugging option */
380 #undef Sudden_Underflow
381 #endif
382 #endif
383
384 #ifndef Flt_Rounds
385 #ifdef FLT_ROUNDS
386 #define Flt_Rounds FLT_ROUNDS
387 #else
388 #define Flt_Rounds 1
389 #endif
390 #endif /*Flt_Rounds*/
391
392 #ifdef Honor_FLT_ROUNDS
393 #undef Check_FLT_ROUNDS
394 #define Check_FLT_ROUNDS
395 #else
396 #define Rounding Flt_Rounds
397 #endif
398
399 #else /* ifndef IEEE_Arith */
400 #undef Check_FLT_ROUNDS
401 #undef Honor_FLT_ROUNDS
402 #undef SET_INEXACT
403 #undef Sudden_Underflow
404 #define Sudden_Underflow
405 #ifdef IBM
406 #undef Flt_Rounds
407 #define Flt_Rounds 0
408 #define Exp_shift 24
409 #define Exp_shift1 24
410 #define Exp_msk1 0x1000000
411 #define Exp_msk11 0x1000000
412 #define Exp_mask 0x7f000000
413 #define P 14
414 #define Nbits 56
415 #define Bias 65
416 #define Emax 248
417 #define Emin (-260)
418 #define Exp_1 0x41000000
419 #define Exp_11 0x41000000
420 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
421 #define Frac_mask 0xffffff
422 #define Frac_mask1 0xffffff
423 #define Bletch 4
424 #define Ten_pmax 22
425 #define Bndry_mask 0xefffff
426 #define Bndry_mask1 0xffffff
427 #define LSB 1
428 #define Sign_bit 0x80000000
429 #define Log2P 4
430 #define Tiny0 0x100000
431 #define Tiny1 0
432 #define Quick_max 14
433 #define Int_max 15
434 #else /* VAX */
435 #undef Flt_Rounds
436 #define Flt_Rounds 1
437 #define Exp_shift 23
438 #define Exp_shift1 7
439 #define Exp_msk1 0x80
440 #define Exp_msk11 0x800000
441 #define Exp_mask 0x7f80
442 #define P 56
443 #define Nbits 56
444 #define Bias 129
445 #define Emax 126
446 #define Emin (-129)
447 #define Exp_1 0x40800000
448 #define Exp_11 0x4080
449 #define Ebits 8
450 #define Frac_mask 0x7fffff
451 #define Frac_mask1 0xffff007f
452 #define Ten_pmax 24
453 #define Bletch 2
454 #define Bndry_mask 0xffff007f
455 #define Bndry_mask1 0xffff007f
456 #define LSB 0x10000
457 #define Sign_bit 0x8000
458 #define Log2P 1
459 #define Tiny0 0x80
460 #define Tiny1 0
461 #define Quick_max 15
462 #define Int_max 15
463 #endif /* IBM, VAX */
464 #endif /* IEEE_Arith */
465
466 #ifndef IEEE_Arith
467 #define ROUND_BIASED
468 #else
469 #ifdef ROUND_BIASED_without_Round_Up
470 #undef ROUND_BIASED
471 #define ROUND_BIASED
472 #endif
473 #endif
474
475 #ifdef RND_PRODQUOT
476 #define rounded_product(a,b) a = rnd_prod(a, b)
477 #define rounded_quotient(a,b) a = rnd_quot(a, b)
478 #ifdef KR_headers
479 extern double rnd_prod(), rnd_quot();
480 #else
481 extern double rnd_prod(double, double), rnd_quot(double, double);
482 #endif
483 #else
484 #define rounded_product(a,b) a *= b
485 #define rounded_quotient(a,b) a /= b
486 #endif
487
488 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
489 #define Big1 0xffffffff
490
491 #ifndef Pack_32
492 #define Pack_32
493 #endif
494
495 typedef struct BCinfo BCinfo;
496 struct
497 BCinfo { int dp0, dp1, dplen, dsign, e0, inexact, nd, nd0, rounding, scale, uflchk; };
498
499 #ifdef KR_headers
500 #define FFFFFFFF ((((unsigned long)0xffff)<<16)|(unsigned long)0xffff)
501 #else
502 #define FFFFFFFF 0xffffffffUL
503 #endif
504
505 #ifdef NO_LONG_LONG
506 #undef ULLong
507 #ifdef Just_16
508 #undef Pack_32
509 /* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
510 * This makes some inner loops simpler and sometimes saves work
511 * during multiplications, but it often seems to make things slightly
512 * slower. Hence the default is now to store 32 bits per Long.
513 */
514 #endif
515 #else /* long long available */
516 #ifndef Llong
517 #define Llong long long
518 #endif
519 #ifndef ULLong
520 #define ULLong unsigned Llong
521 #endif
522 #endif /* NO_LONG_LONG */
523
524 #ifndef MULTIPLE_THREADS
525 #define ACQUIRE_DTOA_LOCK(n) /*nothing*/
526 #define FREE_DTOA_LOCK(n) /*nothing*/
527 #endif
528
529 #define Kmax 7
530
531 #ifdef __cplusplus
532 extern "C" double strtod(const char *s00, char **se);
533 extern "C" char *dtoa(double d, int mode, int ndigits,
534 int *decpt, int *sign, char **rve);
535 #endif
536
537 struct
538 Bigint {
539 struct Bigint *next;
540 int k, maxwds, sign, wds;
541 ULong x[1];
542 };
543
544 typedef struct Bigint Bigint;
545
546 static Bigint *freelist[Kmax+1];
547
548 static void destroy_freelist(void);
549
550 #ifdef ZTS
551 static MUTEX_T dtoa_mutex;
552 static MUTEX_T pow5mult_mutex;
553 #endif /* ZTS */
554
zend_startup_strtod(void)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 /* }}} */
zend_shutdown_strtod(void)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 *
Balloc(k)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
Bfree(v)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 *
multadd(b,m,a)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 *
s2b(s,nd0,nd,y9,dplen)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
hi0bits(x)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
lo0bits(y)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 *
i2b(i)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 *
mult(a,b)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 *
pow5mult(b,k)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 *
lshift(b,k)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
cmp(a,b)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 *
diff(a,b)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
ulp(x)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
b2d(a,e)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 *
d2b(d,e,bits)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
ratio(a,b)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
match(sp,t)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
hexnan(rvp,sp)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
increment(b)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
rshift(b,k)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
any_on(b,k)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
gethex(sp,rvp,rounding,sign)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
dshift(b,p2)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
quorem(b,S)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
sulp(x,bc)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
bigcomp(rv,s0,bc)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
zend_strtod(s00,se)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 if (nd < 0) {
2681 /* overflow */
2682 nd = DBL_DIG + 2;
2683 }
2684 if (nf < 0) {
2685 /* overflow */
2686 nf = DBL_DIG + 2;
2687 }
2688 e = 0;
2689 if (c == 'e' || c == 'E') {
2690 if (!nd && !nz && !nz0) {
2691 goto ret0;
2692 }
2693 s00 = s;
2694 esign = 0;
2695 switch(c = *++s) {
2696 case '-':
2697 esign = 1;
2698 case '+':
2699 c = *++s;
2700 }
2701 if (c >= '0' && c <= '9') {
2702 while(c == '0')
2703 c = *++s;
2704 if (c > '0' && c <= '9') {
2705 L = c - '0';
2706 s1 = s;
2707 while((c = *++s) >= '0' && c <= '9')
2708 L = 10*L + c - '0';
2709 if (s - s1 > 8 || L > 19999)
2710 /* Avoid confusion from exponents
2711 * so large that e might overflow.
2712 */
2713 e = 19999; /* safe for 16 bit ints */
2714 else
2715 e = (int)L;
2716 if (esign)
2717 e = -e;
2718 }
2719 else
2720 e = 0;
2721 }
2722 else
2723 s = s00;
2724 }
2725 if (!nd) {
2726 if (!nz && !nz0) {
2727 #ifdef INFNAN_CHECK
2728 /* Check for Nan and Infinity */
2729 if (!bc.dplen)
2730 switch(c) {
2731 case 'i':
2732 case 'I':
2733 if (match(&s,"nf")) {
2734 --s;
2735 if (!match(&s,"inity"))
2736 ++s;
2737 word0(&rv) = 0x7ff00000;
2738 word1(&rv) = 0;
2739 goto ret;
2740 }
2741 break;
2742 case 'n':
2743 case 'N':
2744 if (match(&s, "an")) {
2745 word0(&rv) = NAN_WORD0;
2746 word1(&rv) = NAN_WORD1;
2747 #ifndef No_Hex_NaN
2748 if (*s == '(') /*)*/
2749 hexnan(&rv, &s);
2750 #endif
2751 goto ret;
2752 }
2753 }
2754 #endif /* INFNAN_CHECK */
2755 ret0:
2756 s = s00;
2757 sign = 0;
2758 }
2759 goto ret;
2760 }
2761 bc.e0 = e1 = e -= nf;
2762
2763 /* Now we have nd0 digits, starting at s0, followed by a
2764 * decimal point, followed by nd-nd0 digits. The number we're
2765 * after is the integer represented by those digits times
2766 * 10**e */
2767
2768 if (!nd0)
2769 nd0 = nd;
2770 k = nd < DBL_DIG + 2 ? nd : DBL_DIG + 2;
2771 dval(&rv) = y;
2772 if (k > 9) {
2773 #ifdef SET_INEXACT
2774 if (k > DBL_DIG)
2775 oldinexact = get_inexact();
2776 #endif
2777 dval(&rv) = tens[k - 9] * dval(&rv) + z;
2778 }
2779 bd0 = 0;
2780 if (nd <= DBL_DIG
2781 #ifndef RND_PRODQUOT
2782 #ifndef Honor_FLT_ROUNDS
2783 && Flt_Rounds == 1
2784 #endif
2785 #endif
2786 ) {
2787 if (!e)
2788 goto ret;
2789 #ifndef ROUND_BIASED_without_Round_Up
2790 if (e > 0) {
2791 if (e <= Ten_pmax) {
2792 #ifdef VAX
2793 goto vax_ovfl_check;
2794 #else
2795 #ifdef Honor_FLT_ROUNDS
2796 /* round correctly FLT_ROUNDS = 2 or 3 */
2797 if (sign) {
2798 rv.d = -rv.d;
2799 sign = 0;
2800 }
2801 #endif
2802 /* rv = */ rounded_product(dval(&rv), tens[e]);
2803 goto ret;
2804 #endif
2805 }
2806 i = DBL_DIG - nd;
2807 if (e <= Ten_pmax + i) {
2808 /* A fancier test would sometimes let us do
2809 * this for larger i values.
2810 */
2811 #ifdef Honor_FLT_ROUNDS
2812 /* round correctly FLT_ROUNDS = 2 or 3 */
2813 if (sign) {
2814 rv.d = -rv.d;
2815 sign = 0;
2816 }
2817 #endif
2818 e -= i;
2819 dval(&rv) *= tens[i];
2820 #ifdef VAX
2821 /* VAX exponent range is so narrow we must
2822 * worry about overflow here...
2823 */
2824 vax_ovfl_check:
2825 word0(&rv) -= P*Exp_msk1;
2826 /* rv = */ rounded_product(dval(&rv), tens[e]);
2827 if ((word0(&rv) & Exp_mask)
2828 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
2829 goto ovfl;
2830 word0(&rv) += P*Exp_msk1;
2831 #else
2832 /* rv = */ rounded_product(dval(&rv), tens[e]);
2833 #endif
2834 goto ret;
2835 }
2836 }
2837 #ifndef Inaccurate_Divide
2838 else if (e >= -Ten_pmax) {
2839 #ifdef Honor_FLT_ROUNDS
2840 /* round correctly FLT_ROUNDS = 2 or 3 */
2841 if (sign) {
2842 rv.d = -rv.d;
2843 sign = 0;
2844 }
2845 #endif
2846 /* rv = */ rounded_quotient(dval(&rv), tens[-e]);
2847 goto ret;
2848 }
2849 #endif
2850 #endif /* ROUND_BIASED_without_Round_Up */
2851 }
2852 e1 += nd - k;
2853
2854 #ifdef IEEE_Arith
2855 #ifdef SET_INEXACT
2856 bc.inexact = 1;
2857 if (k <= DBL_DIG)
2858 oldinexact = get_inexact();
2859 #endif
2860 #ifdef Avoid_Underflow
2861 bc.scale = 0;
2862 #endif
2863 #ifdef Honor_FLT_ROUNDS
2864 if (bc.rounding >= 2) {
2865 if (sign)
2866 bc.rounding = bc.rounding == 2 ? 0 : 2;
2867 else
2868 if (bc.rounding != 2)
2869 bc.rounding = 0;
2870 }
2871 #endif
2872 #endif /*IEEE_Arith*/
2873
2874 /* Get starting approximation = rv * 10**e1 */
2875
2876 if (e1 > 0) {
2877 if ((i = e1 & 15))
2878 dval(&rv) *= tens[i];
2879 if (e1 &= ~15) {
2880 if (e1 > DBL_MAX_10_EXP) {
2881 ovfl:
2882 /* Can't trust HUGE_VAL */
2883 #ifdef IEEE_Arith
2884 #ifdef Honor_FLT_ROUNDS
2885 switch(bc.rounding) {
2886 case 0: /* toward 0 */
2887 case 3: /* toward -infinity */
2888 word0(&rv) = Big0;
2889 word1(&rv) = Big1;
2890 break;
2891 default:
2892 word0(&rv) = Exp_mask;
2893 word1(&rv) = 0;
2894 }
2895 #else /*Honor_FLT_ROUNDS*/
2896 word0(&rv) = Exp_mask;
2897 word1(&rv) = 0;
2898 #endif /*Honor_FLT_ROUNDS*/
2899 #ifdef SET_INEXACT
2900 /* set overflow bit */
2901 dval(&rv0) = 1e300;
2902 dval(&rv0) *= dval(&rv0);
2903 #endif
2904 #else /*IEEE_Arith*/
2905 word0(&rv) = Big0;
2906 word1(&rv) = Big1;
2907 #endif /*IEEE_Arith*/
2908 range_err:
2909 if (bd0) {
2910 Bfree(bb);
2911 Bfree(bd);
2912 Bfree(bs);
2913 Bfree(bd0);
2914 Bfree(delta);
2915 }
2916 #ifndef NO_ERRNO
2917 errno = ERANGE;
2918 #endif
2919 goto ret;
2920 }
2921 e1 >>= 4;
2922 for(j = 0; e1 > 1; j++, e1 >>= 1)
2923 if (e1 & 1)
2924 dval(&rv) *= bigtens[j];
2925 /* The last multiplication could overflow. */
2926 word0(&rv) -= P*Exp_msk1;
2927 dval(&rv) *= bigtens[j];
2928 if ((z = word0(&rv) & Exp_mask)
2929 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
2930 goto ovfl;
2931 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
2932 /* set to largest number */
2933 /* (Can't trust DBL_MAX) */
2934 word0(&rv) = Big0;
2935 word1(&rv) = Big1;
2936 }
2937 else
2938 word0(&rv) += P*Exp_msk1;
2939 }
2940 }
2941 else if (e1 < 0) {
2942 e1 = -e1;
2943 if ((i = e1 & 15))
2944 dval(&rv) /= tens[i];
2945 if (e1 >>= 4) {
2946 if (e1 >= 1 << n_bigtens)
2947 goto undfl;
2948 #ifdef Avoid_Underflow
2949 if (e1 & Scale_Bit)
2950 bc.scale = 2*P;
2951 for(j = 0; e1 > 0; j++, e1 >>= 1)
2952 if (e1 & 1)
2953 dval(&rv) *= tinytens[j];
2954 if (bc.scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask)
2955 >> Exp_shift)) > 0) {
2956 /* scaled rv is denormal; clear j low bits */
2957 if (j >= 32) {
2958 if (j > 54)
2959 goto undfl;
2960 word1(&rv) = 0;
2961 if (j >= 53)
2962 word0(&rv) = (P+2)*Exp_msk1;
2963 else
2964 word0(&rv) &= 0xffffffff << (j-32);
2965 }
2966 else
2967 word1(&rv) &= 0xffffffff << j;
2968 }
2969 #else
2970 for(j = 0; e1 > 1; j++, e1 >>= 1)
2971 if (e1 & 1)
2972 dval(&rv) *= tinytens[j];
2973 /* The last multiplication could underflow. */
2974 dval(&rv0) = dval(&rv);
2975 dval(&rv) *= tinytens[j];
2976 if (!dval(&rv)) {
2977 dval(&rv) = 2.*dval(&rv0);
2978 dval(&rv) *= tinytens[j];
2979 #endif
2980 if (!dval(&rv)) {
2981 undfl:
2982 dval(&rv) = 0.;
2983 goto range_err;
2984 }
2985 #ifndef Avoid_Underflow
2986 word0(&rv) = Tiny0;
2987 word1(&rv) = Tiny1;
2988 /* The refinement below will clean
2989 * this approximation up.
2990 */
2991 }
2992 #endif
2993 }
2994 }
2995
2996 /* Now the hard part -- adjusting rv to the correct value.*/
2997
2998 /* Put digits into bd: true value = bd * 10^e */
2999
3000 bc.nd = nd - nz1;
3001 #ifndef NO_STRTOD_BIGCOMP
3002 bc.nd0 = nd0; /* Only needed if nd > strtod_diglim, but done here */
3003 /* to silence an erroneous warning about bc.nd0 */
3004 /* possibly not being initialized. */
3005 if (nd > strtod_diglim) {
3006 /* ASSERT(strtod_diglim >= 18); 18 == one more than the */
3007 /* minimum number of decimal digits to distinguish double values */
3008 /* in IEEE arithmetic. */
3009 i = j = 18;
3010 if (i > nd0)
3011 j += bc.dplen;
3012 for(;;) {
3013 if (--j < bc.dp1 && j >= bc.dp0)
3014 j = bc.dp0 - 1;
3015 if (s0[j] != '0')
3016 break;
3017 --i;
3018 }
3019 e += nd - i;
3020 nd = i;
3021 if (nd0 > nd)
3022 nd0 = nd;
3023 if (nd < 9) { /* must recompute y */
3024 y = 0;
3025 for(i = 0; i < nd0; ++i)
3026 y = 10*y + s0[i] - '0';
3027 for(j = bc.dp1; i < nd; ++i)
3028 y = 10*y + s0[j++] - '0';
3029 }
3030 }
3031 #endif
3032 bd0 = s2b(s0, nd0, nd, y, bc.dplen);
3033
3034 for(;;) {
3035 bd = Balloc(bd0->k);
3036 Bcopy(bd, bd0);
3037 bb = d2b(&rv, &bbe, &bbbits); /* rv = bb * 2^bbe */
3038 bs = i2b(1);
3039
3040 if (e >= 0) {
3041 bb2 = bb5 = 0;
3042 bd2 = bd5 = e;
3043 }
3044 else {
3045 bb2 = bb5 = -e;
3046 bd2 = bd5 = 0;
3047 }
3048 if (bbe >= 0)
3049 bb2 += bbe;
3050 else
3051 bd2 -= bbe;
3052 bs2 = bb2;
3053 #ifdef Honor_FLT_ROUNDS
3054 if (bc.rounding != 1)
3055 bs2++;
3056 #endif
3057 #ifdef Avoid_Underflow
3058 Lsb = LSB;
3059 Lsb1 = 0;
3060 j = bbe - bc.scale;
3061 i = j + bbbits - 1; /* logb(rv) */
3062 j = P + 1 - bbbits;
3063 if (i < Emin) { /* denormal */
3064 i = Emin - i;
3065 j -= i;
3066 if (i < 32)
3067 Lsb <<= i;
3068 else if (i < 52)
3069 Lsb1 = Lsb << (i-32);
3070 else
3071 Lsb1 = Exp_mask;
3072 }
3073 #else /*Avoid_Underflow*/
3074 #ifdef Sudden_Underflow
3075 #ifdef IBM
3076 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
3077 #else
3078 j = P + 1 - bbbits;
3079 #endif
3080 #else /*Sudden_Underflow*/
3081 j = bbe;
3082 i = j + bbbits - 1; /* logb(rv) */
3083 if (i < Emin) /* denormal */
3084 j += P - Emin;
3085 else
3086 j = P + 1 - bbbits;
3087 #endif /*Sudden_Underflow*/
3088 #endif /*Avoid_Underflow*/
3089 bb2 += j;
3090 bd2 += j;
3091 #ifdef Avoid_Underflow
3092 bd2 += bc.scale;
3093 #endif
3094 i = bb2 < bd2 ? bb2 : bd2;
3095 if (i > bs2)
3096 i = bs2;
3097 if (i > 0) {
3098 bb2 -= i;
3099 bd2 -= i;
3100 bs2 -= i;
3101 }
3102 if (bb5 > 0) {
3103 bs = pow5mult(bs, bb5);
3104 bb1 = mult(bs, bb);
3105 Bfree(bb);
3106 bb = bb1;
3107 }
3108 if (bb2 > 0)
3109 bb = lshift(bb, bb2);
3110 if (bd5 > 0)
3111 bd = pow5mult(bd, bd5);
3112 if (bd2 > 0)
3113 bd = lshift(bd, bd2);
3114 if (bs2 > 0)
3115 bs = lshift(bs, bs2);
3116 delta = diff(bb, bd);
3117 bc.dsign = delta->sign;
3118 delta->sign = 0;
3119 i = cmp(delta, bs);
3120 #ifndef NO_STRTOD_BIGCOMP /*{*/
3121 if (bc.nd > nd && i <= 0) {
3122 if (bc.dsign) {
3123 /* Must use bigcomp(). */
3124 req_bigcomp = 1;
3125 break;
3126 }
3127 #ifdef Honor_FLT_ROUNDS
3128 if (bc.rounding != 1) {
3129 if (i < 0) {
3130 req_bigcomp = 1;
3131 break;
3132 }
3133 }
3134 else
3135 #endif
3136 i = -1; /* Discarded digits make delta smaller. */
3137 }
3138 #endif /*}*/
3139 #ifdef Honor_FLT_ROUNDS /*{*/
3140 if (bc.rounding != 1) {
3141 if (i < 0) {
3142 /* Error is less than an ulp */
3143 if (!delta->x[0] && delta->wds <= 1) {
3144 /* exact */
3145 #ifdef SET_INEXACT
3146 bc.inexact = 0;
3147 #endif
3148 break;
3149 }
3150 if (bc.rounding) {
3151 if (bc.dsign) {
3152 adj.d = 1.;
3153 goto apply_adj;
3154 }
3155 }
3156 else if (!bc.dsign) {
3157 adj.d = -1.;
3158 if (!word1(&rv)
3159 && !(word0(&rv) & Frac_mask)) {
3160 y = word0(&rv) & Exp_mask;
3161 #ifdef Avoid_Underflow
3162 if (!bc.scale || y > 2*P*Exp_msk1)
3163 #else
3164 if (y)
3165 #endif
3166 {
3167 delta = lshift(delta,Log2P);
3168 if (cmp(delta, bs) <= 0)
3169 adj.d = -0.5;
3170 }
3171 }
3172 apply_adj:
3173 #ifdef Avoid_Underflow /*{*/
3174 if (bc.scale && (y = word0(&rv) & Exp_mask)
3175 <= 2*P*Exp_msk1)
3176 word0(&adj) += (2*P+1)*Exp_msk1 - y;
3177 #else
3178 #ifdef Sudden_Underflow
3179 if ((word0(&rv) & Exp_mask) <=
3180 P*Exp_msk1) {
3181 word0(&rv) += P*Exp_msk1;
3182 dval(&rv) += adj.d*ulp(dval(&rv));
3183 word0(&rv) -= P*Exp_msk1;
3184 }
3185 else
3186 #endif /*Sudden_Underflow*/
3187 #endif /*Avoid_Underflow}*/
3188 dval(&rv) += adj.d*ulp(&rv);
3189 }
3190 break;
3191 }
3192 adj.d = ratio(delta, bs);
3193 if (adj.d < 1.)
3194 adj.d = 1.;
3195 if (adj.d <= 0x7ffffffe) {
3196 /* adj = rounding ? ceil(adj) : floor(adj); */
3197 y = adj.d;
3198 if (y != adj.d) {
3199 if (!((bc.rounding>>1) ^ bc.dsign))
3200 y++;
3201 adj.d = y;
3202 }
3203 }
3204 #ifdef Avoid_Underflow /*{*/
3205 if (bc.scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
3206 word0(&adj) += (2*P+1)*Exp_msk1 - y;
3207 #else
3208 #ifdef Sudden_Underflow
3209 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
3210 word0(&rv) += P*Exp_msk1;
3211 adj.d *= ulp(dval(&rv));
3212 if (bc.dsign)
3213 dval(&rv) += adj.d;
3214 else
3215 dval(&rv) -= adj.d;
3216 word0(&rv) -= P*Exp_msk1;
3217 goto cont;
3218 }
3219 #endif /*Sudden_Underflow*/
3220 #endif /*Avoid_Underflow}*/
3221 adj.d *= ulp(&rv);
3222 if (bc.dsign) {
3223 if (word0(&rv) == Big0 && word1(&rv) == Big1)
3224 goto ovfl;
3225 dval(&rv) += adj.d;
3226 }
3227 else
3228 dval(&rv) -= adj.d;
3229 goto cont;
3230 }
3231 #endif /*}Honor_FLT_ROUNDS*/
3232
3233 if (i < 0) {
3234 /* Error is less than half an ulp -- check for
3235 * special case of mantissa a power of two.
3236 */
3237 if (bc.dsign || word1(&rv) || word0(&rv) & Bndry_mask
3238 #ifdef IEEE_Arith /*{*/
3239 #ifdef Avoid_Underflow
3240 || (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1
3241 #else
3242 || (word0(&rv) & Exp_mask) <= Exp_msk1
3243 #endif
3244 #endif /*}*/
3245 ) {
3246 #ifdef SET_INEXACT
3247 if (!delta->x[0] && delta->wds <= 1)
3248 bc.inexact = 0;
3249 #endif
3250 break;
3251 }
3252 if (!delta->x[0] && delta->wds <= 1) {
3253 /* exact result */
3254 #ifdef SET_INEXACT
3255 bc.inexact = 0;
3256 #endif
3257 break;
3258 }
3259 delta = lshift(delta,Log2P);
3260 if (cmp(delta, bs) > 0)
3261 goto drop_down;
3262 break;
3263 }
3264 if (i == 0) {
3265 /* exactly half-way between */
3266 if (bc.dsign) {
3267 if ((word0(&rv) & Bndry_mask1) == Bndry_mask1
3268 && word1(&rv) == (
3269 #ifdef Avoid_Underflow
3270 (bc.scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
3271 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
3272 #endif
3273 0xffffffff)) {
3274 /*boundary case -- increment exponent*/
3275 if (word0(&rv) == Big0 && word1(&rv) == Big1)
3276 goto ovfl;
3277 word0(&rv) = (word0(&rv) & Exp_mask)
3278 + Exp_msk1
3279 #ifdef IBM
3280 | Exp_msk1 >> 4
3281 #endif
3282 ;
3283 word1(&rv) = 0;
3284 #ifdef Avoid_Underflow
3285 bc.dsign = 0;
3286 #endif
3287 break;
3288 }
3289 }
3290 else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) {
3291 drop_down:
3292 /* boundary case -- decrement exponent */
3293 #ifdef Sudden_Underflow /*{{*/
3294 L = word0(&rv) & Exp_mask;
3295 #ifdef IBM
3296 if (L < Exp_msk1)
3297 #else
3298 #ifdef Avoid_Underflow
3299 if (L <= (bc.scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
3300 #else
3301 if (L <= Exp_msk1)
3302 #endif /*Avoid_Underflow*/
3303 #endif /*IBM*/
3304 {
3305 if (bc.nd >nd) {
3306 bc.uflchk = 1;
3307 break;
3308 }
3309 goto undfl;
3310 }
3311 L -= Exp_msk1;
3312 #else /*Sudden_Underflow}{*/
3313 #ifdef Avoid_Underflow
3314 if (bc.scale) {
3315 L = word0(&rv) & Exp_mask;
3316 if (L <= (2*P+1)*Exp_msk1) {
3317 if (L > (P+2)*Exp_msk1)
3318 /* round even ==> */
3319 /* accept rv */
3320 break;
3321 /* rv = smallest denormal */
3322 if (bc.nd >nd) {
3323 bc.uflchk = 1;
3324 break;
3325 }
3326 goto undfl;
3327 }
3328 }
3329 #endif /*Avoid_Underflow*/
3330 L = (word0(&rv) & Exp_mask) - Exp_msk1;
3331 #endif /*Sudden_Underflow}}*/
3332 word0(&rv) = L | Bndry_mask1;
3333 word1(&rv) = 0xffffffff;
3334 #ifdef IBM
3335 goto cont;
3336 #else
3337 #ifndef NO_STRTOD_BIGCOMP
3338 if (bc.nd > nd)
3339 goto cont;
3340 #endif
3341 break;
3342 #endif
3343 }
3344 #ifndef ROUND_BIASED
3345 #ifdef Avoid_Underflow
3346 if (Lsb1) {
3347 if (!(word0(&rv) & Lsb1))
3348 break;
3349 }
3350 else if (!(word1(&rv) & Lsb))
3351 break;
3352 #else
3353 if (!(word1(&rv) & LSB))
3354 break;
3355 #endif
3356 #endif
3357 if (bc.dsign)
3358 #ifdef Avoid_Underflow
3359 dval(&rv) += sulp(&rv, &bc);
3360 #else
3361 dval(&rv) += ulp(&rv);
3362 #endif
3363 #ifndef ROUND_BIASED
3364 else {
3365 #ifdef Avoid_Underflow
3366 dval(&rv) -= sulp(&rv, &bc);
3367 #else
3368 dval(&rv) -= ulp(&rv);
3369 #endif
3370 #ifndef Sudden_Underflow
3371 if (!dval(&rv)) {
3372 if (bc.nd >nd) {
3373 bc.uflchk = 1;
3374 break;
3375 }
3376 goto undfl;
3377 }
3378 #endif
3379 }
3380 #ifdef Avoid_Underflow
3381 bc.dsign = 1 - bc.dsign;
3382 #endif
3383 #endif
3384 break;
3385 }
3386 if ((aadj = ratio(delta, bs)) <= 2.) {
3387 if (bc.dsign)
3388 aadj = aadj1 = 1.;
3389 else if (word1(&rv) || word0(&rv) & Bndry_mask) {
3390 #ifndef Sudden_Underflow
3391 if (word1(&rv) == Tiny1 && !word0(&rv)) {
3392 if (bc.nd >nd) {
3393 bc.uflchk = 1;
3394 break;
3395 }
3396 goto undfl;
3397 }
3398 #endif
3399 aadj = 1.;
3400 aadj1 = -1.;
3401 }
3402 else {
3403 /* special case -- power of FLT_RADIX to be */
3404 /* rounded down... */
3405
3406 if (aadj < 2./FLT_RADIX)
3407 aadj = 1./FLT_RADIX;
3408 else
3409 aadj *= 0.5;
3410 aadj1 = -aadj;
3411 }
3412 }
3413 else {
3414 aadj *= 0.5;
3415 aadj1 = bc.dsign ? aadj : -aadj;
3416 #ifdef Check_FLT_ROUNDS
3417 switch(bc.rounding) {
3418 case 2: /* towards +infinity */
3419 aadj1 -= 0.5;
3420 break;
3421 case 0: /* towards 0 */
3422 case 3: /* towards -infinity */
3423 aadj1 += 0.5;
3424 }
3425 #else
3426 if (Flt_Rounds == 0)
3427 aadj1 += 0.5;
3428 #endif /*Check_FLT_ROUNDS*/
3429 }
3430 y = word0(&rv) & Exp_mask;
3431
3432 /* Check for overflow */
3433
3434 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
3435 dval(&rv0) = dval(&rv);
3436 word0(&rv) -= P*Exp_msk1;
3437 adj.d = aadj1 * ulp(&rv);
3438 dval(&rv) += adj.d;
3439 if ((word0(&rv) & Exp_mask) >=
3440 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
3441 if (word0(&rv0) == Big0 && word1(&rv0) == Big1)
3442 goto ovfl;
3443 word0(&rv) = Big0;
3444 word1(&rv) = Big1;
3445 goto cont;
3446 }
3447 else
3448 word0(&rv) += P*Exp_msk1;
3449 }
3450 else {
3451 #ifdef Avoid_Underflow
3452 if (bc.scale && y <= 2*P*Exp_msk1) {
3453 if (aadj <= 0x7fffffff) {
3454 if ((z = aadj) <= 0)
3455 z = 1;
3456 aadj = z;
3457 aadj1 = bc.dsign ? aadj : -aadj;
3458 }
3459 dval(&aadj2) = aadj1;
3460 word0(&aadj2) += (2*P+1)*Exp_msk1 - y;
3461 aadj1 = dval(&aadj2);
3462 adj.d = aadj1 * ulp(&rv);
3463 dval(&rv) += adj.d;
3464 if (rv.d == 0.)
3465 #ifdef NO_STRTOD_BIGCOMP
3466 goto undfl;
3467 #else
3468 {
3469 req_bigcomp = 1;
3470 break;
3471 }
3472 #endif
3473 }
3474 else {
3475 adj.d = aadj1 * ulp(&rv);
3476 dval(&rv) += adj.d;
3477 }
3478 #else
3479 #ifdef Sudden_Underflow
3480 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
3481 dval(&rv0) = dval(&rv);
3482 word0(&rv) += P*Exp_msk1;
3483 adj.d = aadj1 * ulp(&rv);
3484 dval(&rv) += adj.d;
3485 #ifdef IBM
3486 if ((word0(&rv) & Exp_mask) < P*Exp_msk1)
3487 #else
3488 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1)
3489 #endif
3490 {
3491 if (word0(&rv0) == Tiny0
3492 && word1(&rv0) == Tiny1) {
3493 if (bc.nd >nd) {
3494 bc.uflchk = 1;
3495 break;
3496 }
3497 goto undfl;
3498 }
3499 word0(&rv) = Tiny0;
3500 word1(&rv) = Tiny1;
3501 goto cont;
3502 }
3503 else
3504 word0(&rv) -= P*Exp_msk1;
3505 }
3506 else {
3507 adj.d = aadj1 * ulp(&rv);
3508 dval(&rv) += adj.d;
3509 }
3510 #else /*Sudden_Underflow*/
3511 /* Compute adj so that the IEEE rounding rules will
3512 * correctly round rv + adj in some half-way cases.
3513 * If rv * ulp(rv) is denormalized (i.e.,
3514 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
3515 * trouble from bits lost to denormalization;
3516 * example: 1.2e-307 .
3517 */
3518 if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
3519 aadj1 = (double)(int)(aadj + 0.5);
3520 if (!bc.dsign)
3521 aadj1 = -aadj1;
3522 }
3523 adj.d = aadj1 * ulp(&rv);
3524 dval(&rv) += adj.d;
3525 #endif /*Sudden_Underflow*/
3526 #endif /*Avoid_Underflow*/
3527 }
3528 z = word0(&rv) & Exp_mask;
3529 #ifndef SET_INEXACT
3530 if (bc.nd == nd) {
3531 #ifdef Avoid_Underflow
3532 if (!bc.scale)
3533 #endif
3534 if (y == z) {
3535 /* Can we stop now? */
3536 L = (Long)aadj;
3537 aadj -= L;
3538 /* The tolerances below are conservative. */
3539 if (bc.dsign || word1(&rv) || word0(&rv) & Bndry_mask) {
3540 if (aadj < .4999999 || aadj > .5000001)
3541 break;
3542 }
3543 else if (aadj < .4999999/FLT_RADIX)
3544 break;
3545 }
3546 }
3547 #endif
3548 cont:
3549 Bfree(bb);
3550 Bfree(bd);
3551 Bfree(bs);
3552 Bfree(delta);
3553 }
3554 Bfree(bb);
3555 Bfree(bd);
3556 Bfree(bs);
3557 Bfree(bd0);
3558 Bfree(delta);
3559 #ifndef NO_STRTOD_BIGCOMP
3560 if (req_bigcomp) {
3561 bd0 = 0;
3562 bc.e0 += nz1;
3563 bigcomp(&rv, s0, &bc);
3564 y = word0(&rv) & Exp_mask;
3565 if (y == Exp_mask)
3566 goto ovfl;
3567 if (y == 0 && rv.d == 0.)
3568 goto undfl;
3569 }
3570 #endif
3571 #ifdef SET_INEXACT
3572 if (bc.inexact) {
3573 if (!oldinexact) {
3574 word0(&rv0) = Exp_1 + (70 << Exp_shift);
3575 word1(&rv0) = 0;
3576 dval(&rv0) += 1.;
3577 }
3578 }
3579 else if (!oldinexact)
3580 clear_inexact();
3581 #endif
3582 #ifdef Avoid_Underflow
3583 if (bc.scale) {
3584 word0(&rv0) = Exp_1 - 2*P*Exp_msk1;
3585 word1(&rv0) = 0;
3586 dval(&rv) *= dval(&rv0);
3587 #ifndef NO_ERRNO
3588 /* try to avoid the bug of testing an 8087 register value */
3589 #ifdef IEEE_Arith
3590 if (!(word0(&rv) & Exp_mask))
3591 #else
3592 if (word0(&rv) == 0 && word1(&rv) == 0)
3593 #endif
3594 errno = ERANGE;
3595 #endif
3596 }
3597 #endif /* Avoid_Underflow */
3598 #ifdef SET_INEXACT
3599 if (bc.inexact && !(word0(&rv) & Exp_mask)) {
3600 /* set underflow bit */
3601 dval(&rv0) = 1e-300;
3602 dval(&rv0) *= dval(&rv0);
3603 }
3604 #endif
3605 ret:
3606 if (se)
3607 *se = (char *)s;
3608 return sign ? -dval(&rv) : dval(&rv);
3609 }
3610
3611 #ifndef MULTIPLE_THREADS
3612 ZEND_TLS char *dtoa_result;
3613 #endif
3614
3615 static char *
3616 #ifdef KR_headers
rv_alloc(i)3617 rv_alloc(i) int i;
3618 #else
3619 rv_alloc(int i)
3620 #endif
3621 {
3622 int j, k, *r;
3623
3624 j = sizeof(ULong);
3625 for(k = 0;
3626 sizeof(Bigint) - sizeof(ULong) - sizeof(int) + (size_t)j <= (size_t)i;
3627 j <<= 1)
3628 k++;
3629 r = (int*)Balloc(k);
3630 *r = k;
3631 return
3632 #ifndef MULTIPLE_THREADS
3633 dtoa_result =
3634 #endif
3635 (char *)(r+1);
3636 }
3637
3638 static char *
3639 #ifdef KR_headers
nrv_alloc(s,rve,n)3640 nrv_alloc(s, rve, n) char *s, **rve; int n;
3641 #else
3642 nrv_alloc(const char *s, char **rve, int n)
3643 #endif
3644 {
3645 char *rv, *t;
3646
3647 t = rv = rv_alloc(n);
3648 while((*t = *s++)) t++;
3649 if (rve)
3650 *rve = t;
3651 return rv;
3652 }
3653
3654 /* freedtoa(s) must be used to free values s returned by dtoa
3655 * when MULTIPLE_THREADS is #defined. It should be used in all cases,
3656 * but for consistency with earlier versions of dtoa, it is optional
3657 * when MULTIPLE_THREADS is not defined.
3658 */
3659
3660 ZEND_API void
3661 #ifdef KR_headers
zend_freedtoa(s)3662 zend_freedtoa(s) char *s;
3663 #else
3664 zend_freedtoa(char *s)
3665 #endif
3666 {
3667 Bigint *b = (Bigint *)((int *)s - 1);
3668 b->maxwds = 1 << (b->k = *(int*)b);
3669 Bfree(b);
3670 #ifndef MULTIPLE_THREADS
3671 if (s == dtoa_result)
3672 dtoa_result = 0;
3673 #endif
3674 }
3675
3676 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
3677 *
3678 * Inspired by "How to Print Floating-Point Numbers Accurately" by
3679 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
3680 *
3681 * Modifications:
3682 * 1. Rather than iterating, we use a simple numeric overestimate
3683 * to determine k = floor(log10(d)). We scale relevant
3684 * quantities using O(log2(k)) rather than O(k) multiplications.
3685 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
3686 * try to generate digits strictly left to right. Instead, we
3687 * compute with fewer bits and propagate the carry if necessary
3688 * when rounding the final digit up. This is often faster.
3689 * 3. Under the assumption that input will be rounded nearest,
3690 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
3691 * That is, we allow equality in stopping tests when the
3692 * round-nearest rule will give the same floating-point value
3693 * as would satisfaction of the stopping test with strict
3694 * inequality.
3695 * 4. We remove common factors of powers of 2 from relevant
3696 * quantities.
3697 * 5. When converting floating-point integers less than 1e16,
3698 * we use floating-point arithmetic rather than resorting
3699 * to multiple-precision integers.
3700 * 6. When asked to produce fewer than 15 digits, we first try
3701 * to get by with floating-point arithmetic; we resort to
3702 * multiple-precision integer arithmetic only if we cannot
3703 * guarantee that the floating-point calculation has given
3704 * the correctly rounded result. For k requested digits and
3705 * "uniformly" distributed input, the probability is
3706 * something like 10^(k-15) that we must resort to the Long
3707 * calculation.
3708 */
3709
3710 ZEND_API char *
zend_dtoa(dd,mode,ndigits,decpt,sign,rve)3711 zend_dtoa
3712 #ifdef KR_headers
3713 (dd, mode, ndigits, decpt, sign, rve)
3714 double dd; int mode, ndigits, *decpt, *sign; char **rve;
3715 #else
3716 (double dd, int mode, int ndigits, int *decpt, int *sign, char **rve)
3717 #endif
3718 {
3719 /* Arguments ndigits, decpt, sign are similar to those
3720 of ecvt and fcvt; trailing zeros are suppressed from
3721 the returned string. If not null, *rve is set to point
3722 to the end of the return value. If d is +-Infinity or NaN,
3723 then *decpt is set to 9999.
3724
3725 mode:
3726 0 ==> shortest string that yields d when read in
3727 and rounded to nearest.
3728 1 ==> like 0, but with Steele & White stopping rule;
3729 e.g. with IEEE P754 arithmetic , mode 0 gives
3730 1e23 whereas mode 1 gives 9.999999999999999e22.
3731 2 ==> max(1,ndigits) significant digits. This gives a
3732 return value similar to that of ecvt, except
3733 that trailing zeros are suppressed.
3734 3 ==> through ndigits past the decimal point. This
3735 gives a return value similar to that from fcvt,
3736 except that trailing zeros are suppressed, and
3737 ndigits can be negative.
3738 4,5 ==> similar to 2 and 3, respectively, but (in
3739 round-nearest mode) with the tests of mode 0 to
3740 possibly return a shorter string that rounds to d.
3741 With IEEE arithmetic and compilation with
3742 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
3743 as modes 2 and 3 when FLT_ROUNDS != 1.
3744 6-9 ==> Debugging modes similar to mode - 4: don't try
3745 fast floating-point estimate (if applicable).
3746
3747 Values of mode other than 0-9 are treated as mode 0.
3748
3749 Sufficient space is allocated to the return value
3750 to hold the suppressed trailing zeros.
3751 */
3752
3753 int bbits, b2, b5, be, dig, i, ieps, ilim = 0, ilim0, ilim1,
3754 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
3755 spec_case = 0, try_quick;
3756 Long L;
3757 #ifndef Sudden_Underflow
3758 int denorm;
3759 ULong x;
3760 #endif
3761 Bigint *b, *b1, *delta, *mlo, *mhi, *S;
3762 U d2, eps, u;
3763 double ds;
3764 char *s, *s0;
3765 #ifndef No_leftright
3766 #ifdef IEEE_Arith
3767 U eps1;
3768 #endif
3769 #endif
3770 #ifdef SET_INEXACT
3771 int inexact, oldinexact;
3772 #endif
3773 #ifdef Honor_FLT_ROUNDS /*{*/
3774 int Rounding;
3775 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
3776 Rounding = Flt_Rounds;
3777 #else /*}{*/
3778 Rounding = 1;
3779 switch(fegetround()) {
3780 case FE_TOWARDZERO: Rounding = 0; break;
3781 case FE_UPWARD: Rounding = 2; break;
3782 case FE_DOWNWARD: Rounding = 3;
3783 }
3784 #endif /*}}*/
3785 #endif /*}*/
3786
3787 #ifndef MULTIPLE_THREADS
3788 if (dtoa_result) {
3789 zend_freedtoa(dtoa_result);
3790 dtoa_result = 0;
3791 }
3792 #endif
3793
3794 u.d = dd;
3795 if (word0(&u) & Sign_bit) {
3796 /* set sign for everything, including 0's and NaNs */
3797 *sign = 1;
3798 word0(&u) &= ~Sign_bit; /* clear sign bit */
3799 }
3800 else
3801 *sign = 0;
3802
3803 #if defined(IEEE_Arith) + defined(VAX)
3804 #ifdef IEEE_Arith
3805 if ((word0(&u) & Exp_mask) == Exp_mask)
3806 #else
3807 if (word0(&u) == 0x8000)
3808 #endif
3809 {
3810 /* Infinity or NaN */
3811 *decpt = 9999;
3812 #ifdef IEEE_Arith
3813 if (!word1(&u) && !(word0(&u) & 0xfffff))
3814 return nrv_alloc("Infinity", rve, 8);
3815 #endif
3816 return nrv_alloc("NaN", rve, 3);
3817 }
3818 #endif
3819 #ifdef IBM
3820 dval(&u) += 0; /* normalize */
3821 #endif
3822 if (!dval(&u)) {
3823 *decpt = 1;
3824 return nrv_alloc("0", rve, 1);
3825 }
3826
3827 #ifdef SET_INEXACT
3828 try_quick = oldinexact = get_inexact();
3829 inexact = 1;
3830 #endif
3831 #ifdef Honor_FLT_ROUNDS
3832 if (Rounding >= 2) {
3833 if (*sign)
3834 Rounding = Rounding == 2 ? 0 : 2;
3835 else
3836 if (Rounding != 2)
3837 Rounding = 0;
3838 }
3839 #endif
3840
3841 b = d2b(&u, &be, &bbits);
3842 #ifdef Sudden_Underflow
3843 i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
3844 #else
3845 if ((i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1)))) {
3846 #endif
3847 dval(&d2) = dval(&u);
3848 word0(&d2) &= Frac_mask1;
3849 word0(&d2) |= Exp_11;
3850 #ifdef IBM
3851 if (j = 11 - hi0bits(word0(&d2) & Frac_mask))
3852 dval(&d2) /= 1 << j;
3853 #endif
3854
3855 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
3856 * log10(x) = log(x) / log(10)
3857 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
3858 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
3859 *
3860 * This suggests computing an approximation k to log10(d) by
3861 *
3862 * k = (i - Bias)*0.301029995663981
3863 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
3864 *
3865 * We want k to be too large rather than too small.
3866 * The error in the first-order Taylor series approximation
3867 * is in our favor, so we just round up the constant enough
3868 * to compensate for any error in the multiplication of
3869 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
3870 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
3871 * adding 1e-13 to the constant term more than suffices.
3872 * Hence we adjust the constant term to 0.1760912590558.
3873 * (We could get a more accurate k by invoking log10,
3874 * but this is probably not worthwhile.)
3875 */
3876
3877 i -= Bias;
3878 #ifdef IBM
3879 i <<= 2;
3880 i += j;
3881 #endif
3882 #ifndef Sudden_Underflow
3883 denorm = 0;
3884 }
3885 else {
3886 /* d is denormalized */
3887
3888 i = bbits + be + (Bias + (P-1) - 1);
3889 x = i > 32 ? word0(&u) << (64 - i) | word1(&u) >> (i - 32)
3890 : word1(&u) << (32 - i);
3891 dval(&d2) = x;
3892 word0(&d2) -= 31*Exp_msk1; /* adjust exponent */
3893 i -= (Bias + (P-1) - 1) + 1;
3894 denorm = 1;
3895 }
3896 #endif
3897 ds = (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
3898 k = (int)ds;
3899 if (ds < 0. && ds != k)
3900 k--; /* want k = floor(ds) */
3901 k_check = 1;
3902 if (k >= 0 && k <= Ten_pmax) {
3903 if (dval(&u) < tens[k])
3904 k--;
3905 k_check = 0;
3906 }
3907 j = bbits - i - 1;
3908 if (j >= 0) {
3909 b2 = 0;
3910 s2 = j;
3911 }
3912 else {
3913 b2 = -j;
3914 s2 = 0;
3915 }
3916 if (k >= 0) {
3917 b5 = 0;
3918 s5 = k;
3919 s2 += k;
3920 }
3921 else {
3922 b2 -= k;
3923 b5 = -k;
3924 s5 = 0;
3925 }
3926 if (mode < 0 || mode > 9)
3927 mode = 0;
3928
3929 #ifndef SET_INEXACT
3930 #ifdef Check_FLT_ROUNDS
3931 try_quick = Rounding == 1;
3932 #else
3933 try_quick = 1;
3934 #endif
3935 #endif /*SET_INEXACT*/
3936
3937 if (mode > 5) {
3938 mode -= 4;
3939 try_quick = 0;
3940 }
3941 leftright = 1;
3942 ilim = ilim1 = -1; /* Values for cases 0 and 1; done here to */
3943 /* silence erroneous "gcc -Wall" warning. */
3944 switch(mode) {
3945 case 0:
3946 case 1:
3947 i = 18;
3948 ndigits = 0;
3949 break;
3950 case 2:
3951 leftright = 0;
3952 /* no break */
3953 case 4:
3954 if (ndigits <= 0)
3955 ndigits = 1;
3956 ilim = ilim1 = i = ndigits;
3957 break;
3958 case 3:
3959 leftright = 0;
3960 /* no break */
3961 case 5:
3962 i = ndigits + k + 1;
3963 ilim = i;
3964 ilim1 = i - 1;
3965 if (i <= 0)
3966 i = 1;
3967 }
3968 s = s0 = rv_alloc(i);
3969
3970 #ifdef Honor_FLT_ROUNDS
3971 if (mode > 1 && Rounding != 1)
3972 leftright = 0;
3973 #endif
3974
3975 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
3976
3977 /* Try to get by with floating-point arithmetic. */
3978
3979 i = 0;
3980 dval(&d2) = dval(&u);
3981 k0 = k;
3982 ilim0 = ilim;
3983 ieps = 2; /* conservative */
3984 if (k > 0) {
3985 ds = tens[k&0xf];
3986 j = k >> 4;
3987 if (j & Bletch) {
3988 /* prevent overflows */
3989 j &= Bletch - 1;
3990 dval(&u) /= bigtens[n_bigtens-1];
3991 ieps++;
3992 }
3993 for(; j; j >>= 1, i++)
3994 if (j & 1) {
3995 ieps++;
3996 ds *= bigtens[i];
3997 }
3998 dval(&u) /= ds;
3999 }
4000 else if ((j1 = -k)) {
4001 dval(&u) *= tens[j1 & 0xf];
4002 for(j = j1 >> 4; j; j >>= 1, i++)
4003 if (j & 1) {
4004 ieps++;
4005 dval(&u) *= bigtens[i];
4006 }
4007 }
4008 if (k_check && dval(&u) < 1. && ilim > 0) {
4009 if (ilim1 <= 0)
4010 goto fast_failed;
4011 ilim = ilim1;
4012 k--;
4013 dval(&u) *= 10.;
4014 ieps++;
4015 }
4016 dval(&eps) = ieps*dval(&u) + 7.;
4017 word0(&eps) -= (P-1)*Exp_msk1;
4018 if (ilim == 0) {
4019 S = mhi = 0;
4020 dval(&u) -= 5.;
4021 if (dval(&u) > dval(&eps))
4022 goto one_digit;
4023 if (dval(&u) < -dval(&eps))
4024 goto no_digits;
4025 goto fast_failed;
4026 }
4027 #ifndef No_leftright
4028 if (leftright) {
4029 /* Use Steele & White method of only
4030 * generating digits needed.
4031 */
4032 dval(&eps) = 0.5/tens[ilim-1] - dval(&eps);
4033 #ifdef IEEE_Arith
4034 if (k0 < 0 && j1 >= 307) {
4035 eps1.d = 1.01e256; /* 1.01 allows roundoff in the next few lines */
4036 word0(&eps1) -= Exp_msk1 * (Bias+P-1);
4037 dval(&eps1) *= tens[j1 & 0xf];
4038 for(i = 0, j = (j1-256) >> 4; j; j >>= 1, i++)
4039 if (j & 1)
4040 dval(&eps1) *= bigtens[i];
4041 if (eps.d < eps1.d)
4042 eps.d = eps1.d;
4043 }
4044 #endif
4045 for(i = 0;;) {
4046 L = dval(&u);
4047 dval(&u) -= L;
4048 *s++ = '0' + (int)L;
4049 if (1. - dval(&u) < dval(&eps))
4050 goto bump_up;
4051 if (dval(&u) < dval(&eps))
4052 goto ret1;
4053 if (++i >= ilim)
4054 break;
4055 dval(&eps) *= 10.;
4056 dval(&u) *= 10.;
4057 }
4058 }
4059 else {
4060 #endif
4061 /* Generate ilim digits, then fix them up. */
4062 dval(&eps) *= tens[ilim-1];
4063 for(i = 1;; i++, dval(&u) *= 10.) {
4064 L = (Long)(dval(&u));
4065 if (!(dval(&u) -= L))
4066 ilim = i;
4067 *s++ = '0' + (int)L;
4068 if (i == ilim) {
4069 if (dval(&u) > 0.5 + dval(&eps))
4070 goto bump_up;
4071 else if (dval(&u) < 0.5 - dval(&eps)) {
4072 while(*--s == '0');
4073 s++;
4074 goto ret1;
4075 }
4076 break;
4077 }
4078 }
4079 #ifndef No_leftright
4080 }
4081 #endif
4082 fast_failed:
4083 s = s0;
4084 dval(&u) = dval(&d2);
4085 k = k0;
4086 ilim = ilim0;
4087 }
4088
4089 /* Do we have a "small" integer? */
4090
4091 if (be >= 0 && k <= Int_max) {
4092 /* Yes. */
4093 ds = tens[k];
4094 if (ndigits < 0 && ilim <= 0) {
4095 S = mhi = 0;
4096 if (ilim < 0 || dval(&u) <= 5*ds)
4097 goto no_digits;
4098 goto one_digit;
4099 }
4100 for(i = 1;; i++, dval(&u) *= 10.) {
4101 L = (Long)(dval(&u) / ds);
4102 dval(&u) -= L*ds;
4103 #ifdef Check_FLT_ROUNDS
4104 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
4105 if (dval(&u) < 0) {
4106 L--;
4107 dval(&u) += ds;
4108 }
4109 #endif
4110 *s++ = '0' + (int)L;
4111 if (!dval(&u)) {
4112 #ifdef SET_INEXACT
4113 inexact = 0;
4114 #endif
4115 break;
4116 }
4117 if (i == ilim) {
4118 #ifdef Honor_FLT_ROUNDS
4119 if (mode > 1)
4120 switch(Rounding) {
4121 case 0: goto ret1;
4122 case 2: goto bump_up;
4123 }
4124 #endif
4125 dval(&u) += dval(&u);
4126 #ifdef ROUND_BIASED
4127 if (dval(&u) >= ds)
4128 #else
4129 if (dval(&u) > ds || (dval(&u) == ds && L & 1))
4130 #endif
4131 {
4132 bump_up:
4133 while(*--s == '9')
4134 if (s == s0) {
4135 k++;
4136 *s = '0';
4137 break;
4138 }
4139 ++*s++;
4140 }
4141 break;
4142 }
4143 }
4144 goto ret1;
4145 }
4146
4147 m2 = b2;
4148 m5 = b5;
4149 mhi = mlo = 0;
4150 if (leftright) {
4151 i =
4152 #ifndef Sudden_Underflow
4153 denorm ? be + (Bias + (P-1) - 1 + 1) :
4154 #endif
4155 #ifdef IBM
4156 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
4157 #else
4158 1 + P - bbits;
4159 #endif
4160 b2 += i;
4161 s2 += i;
4162 mhi = i2b(1);
4163 }
4164 if (m2 > 0 && s2 > 0) {
4165 i = m2 < s2 ? m2 : s2;
4166 b2 -= i;
4167 m2 -= i;
4168 s2 -= i;
4169 }
4170 if (b5 > 0) {
4171 if (leftright) {
4172 if (m5 > 0) {
4173 mhi = pow5mult(mhi, m5);
4174 b1 = mult(mhi, b);
4175 Bfree(b);
4176 b = b1;
4177 }
4178 if ((j = b5 - m5))
4179 b = pow5mult(b, j);
4180 }
4181 else
4182 b = pow5mult(b, b5);
4183 }
4184 S = i2b(1);
4185 if (s5 > 0)
4186 S = pow5mult(S, s5);
4187
4188 /* Check for special case that d is a normalized power of 2. */
4189
4190 spec_case = 0;
4191 if ((mode < 2 || leftright)
4192 #ifdef Honor_FLT_ROUNDS
4193 && Rounding == 1
4194 #endif
4195 ) {
4196 if (!word1(&u) && !(word0(&u) & Bndry_mask)
4197 #ifndef Sudden_Underflow
4198 && word0(&u) & (Exp_mask & ~Exp_msk1)
4199 #endif
4200 ) {
4201 /* The special case */
4202 b2 += Log2P;
4203 s2 += Log2P;
4204 spec_case = 1;
4205 }
4206 }
4207
4208 /* Arrange for convenient computation of quotients:
4209 * shift left if necessary so divisor has 4 leading 0 bits.
4210 *
4211 * Perhaps we should just compute leading 28 bits of S once
4212 * and for all and pass them and a shift to quorem, so it
4213 * can do shifts and ors to compute the numerator for q.
4214 */
4215 i = dshift(S, s2);
4216 b2 += i;
4217 m2 += i;
4218 s2 += i;
4219 if (b2 > 0)
4220 b = lshift(b, b2);
4221 if (s2 > 0)
4222 S = lshift(S, s2);
4223 if (k_check) {
4224 if (cmp(b,S) < 0) {
4225 k--;
4226 b = multadd(b, 10, 0); /* we botched the k estimate */
4227 if (leftright)
4228 mhi = multadd(mhi, 10, 0);
4229 ilim = ilim1;
4230 }
4231 }
4232 if (ilim <= 0 && (mode == 3 || mode == 5)) {
4233 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
4234 /* no digits, fcvt style */
4235 no_digits:
4236 k = -1 - ndigits;
4237 goto ret;
4238 }
4239 one_digit:
4240 *s++ = '1';
4241 k++;
4242 goto ret;
4243 }
4244 if (leftright) {
4245 if (m2 > 0)
4246 mhi = lshift(mhi, m2);
4247
4248 /* Compute mlo -- check for special case
4249 * that d is a normalized power of 2.
4250 */
4251
4252 mlo = mhi;
4253 if (spec_case) {
4254 mhi = Balloc(mhi->k);
4255 Bcopy(mhi, mlo);
4256 mhi = lshift(mhi, Log2P);
4257 }
4258
4259 for(i = 1;;i++) {
4260 dig = quorem(b,S) + '0';
4261 /* Do we yet have the shortest decimal string
4262 * that will round to d?
4263 */
4264 j = cmp(b, mlo);
4265 delta = diff(S, mhi);
4266 j1 = delta->sign ? 1 : cmp(b, delta);
4267 Bfree(delta);
4268 #ifndef ROUND_BIASED
4269 if (j1 == 0 && mode != 1 && !(word1(&u) & 1)
4270 #ifdef Honor_FLT_ROUNDS
4271 && Rounding >= 1
4272 #endif
4273 ) {
4274 if (dig == '9')
4275 goto round_9_up;
4276 if (j > 0)
4277 dig++;
4278 #ifdef SET_INEXACT
4279 else if (!b->x[0] && b->wds <= 1)
4280 inexact = 0;
4281 #endif
4282 *s++ = dig;
4283 goto ret;
4284 }
4285 #endif
4286 if (j < 0 || (j == 0 && mode != 1
4287 #ifndef ROUND_BIASED
4288 && !(word1(&u) & 1)
4289 #endif
4290 )) {
4291 if (!b->x[0] && b->wds <= 1) {
4292 #ifdef SET_INEXACT
4293 inexact = 0;
4294 #endif
4295 goto accept_dig;
4296 }
4297 #ifdef Honor_FLT_ROUNDS
4298 if (mode > 1)
4299 switch(Rounding) {
4300 case 0: goto accept_dig;
4301 case 2: goto keep_dig;
4302 }
4303 #endif /*Honor_FLT_ROUNDS*/
4304 if (j1 > 0) {
4305 b = lshift(b, 1);
4306 j1 = cmp(b, S);
4307 #ifdef ROUND_BIASED
4308 if (j1 >= 0 /*)*/
4309 #else
4310 if ((j1 > 0 || (j1 == 0 && dig & 1))
4311 #endif
4312 && dig++ == '9')
4313 goto round_9_up;
4314 }
4315 accept_dig:
4316 *s++ = dig;
4317 goto ret;
4318 }
4319 if (j1 > 0) {
4320 #ifdef Honor_FLT_ROUNDS
4321 if (!Rounding)
4322 goto accept_dig;
4323 #endif
4324 if (dig == '9') { /* possible if i == 1 */
4325 round_9_up:
4326 *s++ = '9';
4327 goto roundoff;
4328 }
4329 *s++ = dig + 1;
4330 goto ret;
4331 }
4332 #ifdef Honor_FLT_ROUNDS
4333 keep_dig:
4334 #endif
4335 *s++ = dig;
4336 if (i == ilim)
4337 break;
4338 b = multadd(b, 10, 0);
4339 if (mlo == mhi)
4340 mlo = mhi = multadd(mhi, 10, 0);
4341 else {
4342 mlo = multadd(mlo, 10, 0);
4343 mhi = multadd(mhi, 10, 0);
4344 }
4345 }
4346 }
4347 else
4348 for(i = 1;; i++) {
4349 *s++ = dig = quorem(b,S) + '0';
4350 if (!b->x[0] && b->wds <= 1) {
4351 #ifdef SET_INEXACT
4352 inexact = 0;
4353 #endif
4354 goto ret;
4355 }
4356 if (i >= ilim)
4357 break;
4358 b = multadd(b, 10, 0);
4359 }
4360
4361 /* Round off last digit */
4362
4363 #ifdef Honor_FLT_ROUNDS
4364 switch(Rounding) {
4365 case 0: goto trimzeros;
4366 case 2: goto roundoff;
4367 }
4368 #endif
4369 b = lshift(b, 1);
4370 j = cmp(b, S);
4371 #ifdef ROUND_BIASED
4372 if (j >= 0)
4373 #else
4374 if (j > 0 || (j == 0 && dig & 1))
4375 #endif
4376 {
4377 roundoff:
4378 while(*--s == '9')
4379 if (s == s0) {
4380 k++;
4381 *s++ = '1';
4382 goto ret;
4383 }
4384 ++*s++;
4385 }
4386 else {
4387 #ifdef Honor_FLT_ROUNDS
4388 trimzeros:
4389 #endif
4390 while(*--s == '0');
4391 s++;
4392 }
4393 ret:
4394 Bfree(S);
4395 if (mhi) {
4396 if (mlo && mlo != mhi)
4397 Bfree(mlo);
4398 Bfree(mhi);
4399 }
4400 ret1:
4401 #ifdef SET_INEXACT
4402 if (inexact) {
4403 if (!oldinexact) {
4404 word0(&u) = Exp_1 + (70 << Exp_shift);
4405 word1(&u) = 0;
4406 dval(&u) += 1.;
4407 }
4408 }
4409 else if (!oldinexact)
4410 clear_inexact();
4411 #endif
4412 Bfree(b);
4413 *s = 0;
4414 *decpt = k + 1;
4415 if (rve)
4416 *rve = s;
4417 return s0;
4418 }
4419
zend_hex_strtod(const char * str,const char ** endptr)4420 ZEND_API double zend_hex_strtod(const char *str, const char **endptr)
4421 {
4422 const char *s = str;
4423 char c;
4424 int any = 0;
4425 double value = 0;
4426
4427 if (*s == '0' && (s[1] == 'x' || s[1] == 'X')) {
4428 s += 2;
4429 }
4430
4431 while ((c = *s++)) {
4432 if (c >= '0' && c <= '9') {
4433 c -= '0';
4434 } else if (c >= 'A' && c <= 'F') {
4435 c -= 'A' - 10;
4436 } else if (c >= 'a' && c <= 'f') {
4437 c -= 'a' - 10;
4438 } else {
4439 break;
4440 }
4441
4442 any = 1;
4443 value = value * 16 + c;
4444 }
4445
4446 if (endptr != NULL) {
4447 *endptr = any ? s - 1 : str;
4448 }
4449
4450 return value;
4451 }
4452
zend_oct_strtod(const char * str,const char ** endptr)4453 ZEND_API double zend_oct_strtod(const char *str, const char **endptr)
4454 {
4455 const char *s = str;
4456 char c;
4457 double value = 0;
4458 int any = 0;
4459
4460 if (str[0] == '\0') {
4461 if (endptr != NULL) {
4462 *endptr = str;
4463 }
4464 return 0.0;
4465 }
4466
4467 /* skip leading zero */
4468 s++;
4469
4470 while ((c = *s++)) {
4471 if (c < '0' || c > '7') {
4472 /* break and return the current value if the number is not well-formed
4473 * that's what Linux strtol() does
4474 */
4475 break;
4476 }
4477 value = value * 8 + c - '0';
4478 any = 1;
4479 }
4480
4481 if (endptr != NULL) {
4482 *endptr = any ? s - 1 : str;
4483 }
4484
4485 return value;
4486 }
4487
zend_bin_strtod(const char * str,const char ** endptr)4488 ZEND_API double zend_bin_strtod(const char *str, const char **endptr)
4489 {
4490 const char *s = str;
4491 char c;
4492 double value = 0;
4493 int any = 0;
4494
4495 if ('0' == *s && ('b' == s[1] || 'B' == s[1])) {
4496 s += 2;
4497 }
4498
4499 while ((c = *s++)) {
4500 /*
4501 * Verify the validity of the current character as a base-2 digit. In
4502 * the event that an invalid digit is found, halt the conversion and
4503 * return the portion which has been converted thus far.
4504 */
4505 if ('0' == c || '1' == c)
4506 value = value * 2 + c - '0';
4507 else
4508 break;
4509
4510 any = 1;
4511 }
4512
4513 /*
4514 * As with many strtoX implementations, should the subject sequence be
4515 * empty or not well-formed, no conversion is performed and the original
4516 * value of str is stored in *endptr, provided that endptr is not a null
4517 * pointer.
4518 */
4519 if (NULL != endptr) {
4520 *endptr = (char *)(any ? s - 1 : str);
4521 }
4522
4523 return value;
4524 }
4525
destroy_freelist(void)4526 static void destroy_freelist(void)
4527 {
4528 int i;
4529 Bigint *tmp;
4530
4531 ACQUIRE_DTOA_LOCK(0)
4532 for (i = 0; i <= Kmax; i++) {
4533 Bigint **listp = &freelist[i];
4534 while ((tmp = *listp) != NULL) {
4535 *listp = tmp->next;
4536 free(tmp);
4537 }
4538 freelist[i] = NULL;
4539 }
4540 FREE_DTOA_LOCK(0)
4541 }
4542
4543 #ifdef __cplusplus
4544 }
4545 #endif
4546 /*
4547 * Local variables:
4548 * tab-width: 4
4549 * c-basic-offset: 4
4550 * End:
4551 * vim600: sw=4 ts=4 fdm=marker
4552 * vim<600: sw=4 ts=4
4553 */
4554
