1 /****************************************************************
3 * The author of this software is David M. Gay.
5 * Copyright (c) 1991 by AT&T.
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.
13 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
14 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY
15 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
16 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
18 ***************************************************************/
20 /* Please send bug reports to
22 AT&T Bell Laboratories, Room 2C-463
24 Murray Hill, NJ 07974-2070
26 dmg@research.att.com or research!dmg
29 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
31 * This strtod returns a nearest machine number to the input decimal
32 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are
33 * broken by the IEEE round-even rule. Otherwise ties are broken by
34 * biased rounding (add half and chop).
36 * Inspired loosely by William D. Clinger's paper "How to Read Floating
37 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
41 * 1. We only require IEEE, IBM, or VAX double-precision
42 * arithmetic (not IEEE double-extended).
43 * 2. We get by with floating-point arithmetic in a case that
44 * Clinger missed -- when we're computing d * 10^n
45 * for a small integer d and the integer n is not too
46 * much larger than 22 (the maximum integer k for which
47 * we can represent 10^k exactly), we may be able to
48 * compute (d*10^k) * 10^(e-k) with just one roundoff.
49 * 3. Rather than a bit-at-a-time adjustment of the binary
50 * result in the hard case, we use floating-point
51 * arithmetic to determine the adjustment to within
52 * one bit; only in really hard cases do we need to
53 * compute a second residual.
54 * 4. Because of 3., we don't need a large table of powers of 10
55 * for ten-to-e (just some small tables, e.g. of 10^k
60 * #define IEEE_LITTLE_ENDIAN for IEEE-arithmetic machines where the least
61 * significant byte has the lowest address.
62 * #define IEEE_BIG_ENDIAN for IEEE-arithmetic machines where the most
63 * significant byte has the lowest address.
64 * #define Long int on machines with 32-bit ints and 64-bit longs.
65 * #define Sudden_Underflow for IEEE-format machines without gradual
66 * underflow (i.e., that flush to zero on underflow).
67 * #define IBM for IBM mainframe-style floating-point arithmetic.
68 * #define VAX for VAX-style floating-point arithmetic.
69 * #define Unsigned_Shifts if >> does treats its left operand as unsigned.
70 * #define No_leftright to omit left-right logic in fast floating-point
71 * computation of dtoa.
72 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
73 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
74 * that use extended-precision instructions to compute rounded
75 * products and quotients) with IBM.
76 * #define ROUND_BIASED for IEEE-format with biased rounding.
77 * #define Inaccurate_Divide for IEEE-format with correctly rounded
78 * products but inaccurate quotients, e.g., for Intel i860.
79 * #define Just_16 to store 16 bits per 32-bit Long when doing high-precision
80 * integer arithmetic. Whether this speeds things up or slows things
81 * down depends on the machine and the number being converted.
82 * #define KR_headers for old-style C function headers.
83 * #define Bad_float_h if your system lacks a float.h or if it does not
84 * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
85 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
86 * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
87 * if memory is available and otherwise does something you deem
88 * appropriate. If MALLOC is undefined, malloc will be invoked
89 * directly -- and assumed always to succeed.
92 #include "hphp/zend/zend-strtod.h"
94 #include "hphp/util/assertions.h"
95 #include "hphp/util/rds-local.h"
103 ///////////////////////////////////////////////////////////////////////////////
105 #if (defined(__APPLE__) || defined(__APPLE_CC__)) && (defined(__BIG_ENDIAN__) || defined(__LITTLE_ENDIAN__))
106 # if defined(__LITTLE_ENDIAN__)
107 # undef WORDS_BIGENDIAN
109 # if defined(__BIG_ENDIAN__)
110 # define WORDS_BIGENDIAN
115 #ifdef WORDS_BIGENDIAN
116 #define IEEE_BIG_ENDIAN
118 #define IEEE_LITTLE_ENDIAN
121 #if defined(__arm__) && !defined(__VFP_FP__)
123 * * Although the CPU is little endian the FP has different
124 * * byte and word endianness. The byte order is still little endian
125 * * but the word order is big endian.
127 #define IEEE_BIG_ENDIAN
128 #undef IEEE_LITTLE_ENDIAN
135 #if defined(_MSC_VER)
136 #define int32_t __int32
137 #define uint32_t unsigned __int32
138 #define IEEE_LITTLE_ENDIAN
142 #define ULong uint32_t
146 extern char *MALLOC();
148 extern void *MALLOC(size_t);
151 #define MALLOC malloc
162 #ifdef IEEE_BIG_ENDIAN
163 #define IEEE_ARITHMETIC
165 #ifdef IEEE_LITTLE_ENDIAN
166 #define IEEE_ARITHMETIC
169 #ifdef IEEE_ARITHMETIC
171 #define DBL_MAX_10_EXP 308
172 #define DBL_MAX_EXP 1024
175 #define DBL_MAX 1.7976931348623157e+308
180 #define DBL_MAX_10_EXP 75
181 #define DBL_MAX_EXP 63
184 #define DBL_MAX 7.2370055773322621e+75
189 #define DBL_MAX_10_EXP 38
190 #define DBL_MAX_EXP 127
193 #define DBL_MAX 1.7014118346046923e+38
198 #define LONG_MAX 2147483647
217 #define CONST /* blank */
223 #ifdef Unsigned_Shifts
224 #define Sign_Extend(a,b) if (b < 0) a |= 0xffff0000;
226 #define Sign_Extend(a,b) /*no-op*/
229 #if defined(IEEE_LITTLE_ENDIAN) + defined(IEEE_BIG_ENDIAN) + defined(VAX) + \
231 #error Exactly one of IEEE_LITTLE_ENDIAN, IEEE_BIG_ENDIAN, VAX, or IBM \
239 #define value(x) ((x).d)
240 #ifdef IEEE_LITTLE_ENDIAN
241 #define word0(x) ((x).ul[1])
242 #define word1(x) ((x).ul[0])
244 #define word0(x) ((x).ul[0])
245 #define word1(x) ((x).ul[1])
248 /* The following definition of Storeinc is appropriate for MIPS processors.
249 * An alternative that might be better on some machines is
250 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
252 #if defined(IEEE_LITTLE_ENDIAN) + defined(VAX) + defined(__arm__)
253 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
254 ((unsigned short *)a)[0] = (unsigned short)c, a++)
256 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
257 ((unsigned short *)a)[1] = (unsigned short)c, a++)
260 /* #define P DBL_MANT_DIG */
261 /* Ten_pmax = floor(P*log(2)/log(5)) */
262 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
263 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
264 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
266 #if defined(IEEE_LITTLE_ENDIAN) + defined(IEEE_BIG_ENDIAN)
268 #define Exp_shift1 20
269 #define Exp_msk1 0x100000
270 #define Exp_msk11 0x100000
271 #define Exp_mask 0x7ff00000
276 #define Exp_1 0x3ff00000
277 #define Exp_11 0x3ff00000
279 #define Frac_mask 0xfffff
280 #define Frac_mask1 0xfffff
283 #define Bndry_mask 0xfffff
284 #define Bndry_mask1 0xfffff
286 #define Sign_bit 0x80000000
292 #define Infinite(x) (word0(x) == 0x7ff00000) /* sufficient test for here */
294 #undef Sudden_Underflow
295 #define Sudden_Underflow
298 #define Exp_shift1 24
299 #define Exp_msk1 0x1000000
300 #define Exp_msk11 0x1000000
301 #define Exp_mask 0x7f000000
304 #define Exp_1 0x41000000
305 #define Exp_11 0x41000000
306 #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
307 #define Frac_mask 0xffffff
308 #define Frac_mask1 0xffffff
311 #define Bndry_mask 0xefffff
312 #define Bndry_mask1 0xffffff
314 #define Sign_bit 0x80000000
316 #define Tiny0 0x100000
323 #define Exp_msk1 0x80
324 #define Exp_msk11 0x800000
325 #define Exp_mask 0x7f80
328 #define Exp_1 0x40800000
329 #define Exp_11 0x4080
331 #define Frac_mask 0x7fffff
332 #define Frac_mask1 0xffff007f
335 #define Bndry_mask 0xffff007f
336 #define Bndry_mask1 0xffff007f
338 #define Sign_bit 0x8000
352 #define rounded_product(a,b) a = rnd_prod(a, b)
353 #define rounded_quotient(a,b) a = rnd_quot(a, b)
355 extern double rnd_prod(), rnd_quot();
357 extern double rnd_prod(double, double), rnd_quot(double, double);
360 #define rounded_product(a,b) a *= b
361 #define rounded_quotient(a,b) a /= b
364 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
365 #define Big1 0xffffffff
368 /* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
369 * * This makes some inner loops simpler and sometimes saves work
370 * * during multiplications, but it often seems to make things slightly
371 * * slower. Hence the default is now to store 32 bits per Long.
384 int k
, maxwds
, sign
, wds
;
388 typedef struct Bigint Bigint
;
392 void destroy_freelist(Bigint
** freelist
);
395 BigintData() : p5s(nullptr) {
396 freelist
= (Bigint
**)calloc(Kmax
+ 1, sizeof(Bigint
*));
399 destroy_freelist(freelist
);
406 static RDS_LOCAL_NO_CHECK(BigintData
, s_bigint_data
);
408 // NOTE: If this has not been called, various functions in this file,
409 // and in other files (e.g. "zend-printf.cpp"), will crash when called.
410 void zend_get_bigint_data() {
411 s_bigint_data
.getCheck();
414 static Bigint
* Balloc(int k
)
421 Bigint
**&freelist
= s_bigint_data
->freelist
;
422 if ((rv
= freelist
[k
])) {
423 freelist
[k
] = rv
->next
;
426 rv
= (Bigint
*)MALLOC(sizeof(Bigint
) + (x
-1)*sizeof(Long
));
427 assertx(rv
!= nullptr);
431 rv
->sign
= rv
->wds
= 0;
435 static void Bfree(Bigint
*v
)
438 Bigint
**&freelist
= s_bigint_data
->freelist
;
439 v
->next
= freelist
[v
->k
];
444 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
445 y->wds*sizeof(Long) + 2*sizeof(int))
447 /* return value is only used as a simple string, so mis-aligned parts
448 * inside the Bigint are not at risk on strict align architectures
450 static char * rv_alloc(int i
) {
455 (int)(sizeof(Bigint
) - sizeof(ULong
) - sizeof(int) + j
) <= i
;
461 return (char *)(r
+1);
465 static char * nrv_alloc(const char *s
, char **rve
, int n
)
469 t
= rv
= rv_alloc(n
);
470 while((*t
= *s
++) !=0) {
479 static Bigint
* multadd(Bigint
*b
, int m
, int a
) /* multiply by m and add a */
494 y
= (xi
& 0xffff) * m
+ a
;
495 z
= (xi
>> 16) * m
+ (y
>> 16);
497 *x
++ = (z
<< 16) + (y
& 0xffff);
506 if (wds
>= b
->maxwds
) {
518 static int hi0bits(ULong x
)
522 if (!(x
& 0xffff0000)) {
526 if (!(x
& 0xff000000)) {
530 if (!(x
& 0xf0000000)) {
534 if (!(x
& 0xc0000000)) {
538 if (!(x
& 0x80000000)) {
540 if (!(x
& 0x40000000)) {
547 static int lo0bits(ULong
*y
)
591 static Bigint
* i2b(int i
)
601 static Bigint
* mult(Bigint
*a
, Bigint
*b
)
606 ULong
*x
, *xa
, *xae
, *xb
, *xbe
, *xc
, *xc0
;
611 if (a
->wds
< b
->wds
) {
620 if (wc
> a
->maxwds
) {
624 for(x
= c
->x
, xa
= x
+ wc
; x
< xa
; x
++) {
633 for(; xb
< xbe
; xb
++, xc0
++) {
634 if ((y
= *xb
& 0xffff)) {
639 z
= (*x
& 0xffff) * y
+ (*xc
& 0xffff) + carry
;
641 z2
= (*x
++ >> 16) * y
+ (*xc
>> 16) + carry
;
648 if ((y
= *xb
>> 16)) {
654 z
= (*x
& 0xffff) * y
+ (*xc
>> 16) + carry
;
657 z2
= (*x
++ >> 16) * y
+ (*xc
& 0xffff) + carry
;
665 for(; xb
< xbe
; xc0
++) {
671 z
= *x
++ * y
+ *xc
+ carry
;
680 for(xc0
= c
->x
, xc
= xc0
+ wc
; wc
> 0 && !*--xc
; --wc
) ;
685 static Bigint
* s2b (CONST
char *s
, int nd0
, int nd
, ULong y9
)
692 for(k
= 0, y
= 1; x
> y
; y
<<= 1, k
++) ;
699 b
->x
[0] = y9
& 0xffff;
700 b
->wds
= (b
->x
[1] = y9
>> 16) ? 2 : 1;
706 do b
= multadd(b
, 10, *s
++ - '0');
713 b
= multadd(b
, 10, *s
++ - '0');
718 static Bigint
* pow5mult(Bigint
*b
, int k
)
720 Bigint
*b1
, *p5
, *p51
;
722 static int p05
[3] = { 5, 25, 125 };
724 Bigint
*&p5s
= s_bigint_data
->p5s
;
726 b
= multadd(b
, p05
[i
-1], 0);
746 if (!(p51
= p5
->next
)) {
747 if (!(p51
= p5
->next
)) {
748 p51
= p5
->next
= mult(p5
,p5
);
758 static Bigint
*lshift(Bigint
*b
, int k
)
762 ULong
*x
, *x1
, *xe
, z
;
771 for(i
= b
->maxwds
; n1
> i
; i
<<= 1) {
776 for(i
= 0; i
< n
; i
++) {
799 *x1
++ = *x
<< k
& 0xffff | z
;
816 static int cmp(Bigint
*a
, Bigint
*b
)
818 ULong
*xa
, *xa0
, *xb
, *xb0
;
831 return *xa
< *xb
? -1 : 1;
839 static Bigint
* diff(Bigint
*a
, Bigint
*b
)
843 Long borrow
, y
; /* We need signed shifts here. */
844 ULong
*xa
, *xae
, *xb
, *xbe
, *xc
;
876 y
= (*xa
& 0xffff) - (*xb
& 0xffff) + borrow
;
878 Sign_Extend(borrow
, y
);
879 z
= (*xa
++ >> 16) - (*xb
++ >> 16) + borrow
;
881 Sign_Extend(borrow
, z
);
885 y
= (*xa
& 0xffff) + borrow
;
887 Sign_Extend(borrow
, y
);
888 z
= (*xa
++ >> 16) + borrow
;
890 Sign_Extend(borrow
, z
);
895 y
= *xa
++ - *xb
++ + borrow
;
897 Sign_Extend(borrow
, y
);
903 Sign_Extend(borrow
, y
);
914 static double ulp (double _x
)
921 L
= (word0(x
) & Exp_mask
) - (P
-1)*Exp_msk1
;
922 #ifndef Sudden_Underflow
930 #ifndef Sudden_Underflow
935 word0(a
) = 0x80000 >> L
;
941 word1(a
) = L
>= 31 ? 1 : 1 << (31 - L
);
951 (a
, e
) Bigint
*a
; int *e
;
956 ULong
*xa
, *xa0
, w
, y
, z
;
973 d0
= Exp_1
| y
>> (Ebits
- k
);
974 w
= xa
> xa0
? *--xa
: 0;
975 d1
= y
<< ((32-Ebits
) + k
) | w
>> (Ebits
- k
);
978 z
= xa
> xa0
? *--xa
: 0;
980 d0
= Exp_1
| y
<< k
| z
>> (32 - k
);
981 y
= xa
> xa0
? *--xa
: 0;
982 d1
= z
<< k
| y
>> (32 - k
);
989 if (k
< Ebits
+ 16) {
990 z
= xa
> xa0
? *--xa
: 0;
991 d0
= Exp_1
| y
<< k
- Ebits
| z
>> Ebits
+ 16 - k
;
992 w
= xa
> xa0
? *--xa
: 0;
993 y
= xa
> xa0
? *--xa
: 0;
994 d1
= z
<< k
+ 16 - Ebits
| w
<< k
- Ebits
| y
>> 16 + Ebits
- k
;
997 z
= xa
> xa0
? *--xa
: 0;
998 w
= xa
> xa0
? *--xa
: 0;
1000 d0
= Exp_1
| y
<< k
+ 16 | z
<< k
| w
>> 16 - k
;
1001 y
= xa
> xa0
? *--xa
: 0;
1002 d1
= w
<< k
+ 16 | y
<< k
;
1006 word0(d
) = d0
>> 16 | d0
<< 16;
1007 word1(d
) = d1
>> 16 | d1
<< 16;
1016 static Bigint
* d2b(double _d
, int *e
, int *bits
)
1028 d0
= word0(d
) >> 16 | word0(d
) << 16;
1029 d1
= word1(d
) >> 16 | word1(d
) << 16;
1043 d0
&= 0x7fffffff; /* clear sign bit, which we ignore */
1044 #ifdef Sudden_Underflow
1045 de
= (int)(d0
>> Exp_shift
);
1050 if ((de
= (int)(d0
>> Exp_shift
)))
1055 if ((k
= lo0bits(&y
))) {
1056 x
[0] = y
| (z
<< (32 - k
));
1061 i
= b
->wds
= (x
[1] = z
) ? 2 : 1;
1070 if (k
= lo0bits(&y
)) {
1072 x
[0] = y
| z
<< 32 - k
& 0xffff;
1073 x
[1] = z
>> k
- 16 & 0xffff;
1078 x
[1] = y
>> 16 | z
<< 16 - k
& 0xffff;
1079 x
[2] = z
>> k
& 0xffff;
1106 #ifndef Sudden_Underflow
1110 *e
= (de
- Bias
- (P
-1) << 2) + k
;
1111 *bits
= 4*P
+ 8 - k
- hi0bits(word0(d
) & Frac_mask
);
1113 *e
= de
- Bias
- (P
-1) + k
;
1116 #ifndef Sudden_Underflow
1118 *e
= de
- Bias
- (P
-1) + 1 + k
;
1120 *bits
= 32*i
- hi0bits(x
[i
-1]);
1122 *bits
= (i
+2)*16 - hi0bits(x
[i
]);
1132 static double ratio (Bigint
*a
, Bigint
*b
)
1137 value(da
) = b2d(a
, &ka
);
1138 value(db
) = b2d(b
, &kb
);
1140 k
= ka
- kb
+ 32*(a
->wds
- b
->wds
);
1142 k
= ka
- kb
+ 16*(a
->wds
- b
->wds
);
1146 word0(da
) += (k
>> 2)*Exp_msk1
;
1152 word0(db
) += (k
>> 2)*Exp_msk1
;
1158 word0(da
) += k
*Exp_msk1
;
1161 word0(db
) += k
*Exp_msk1
;
1164 return value(da
) / value(db
);
1169 1e0
, 1e1
, 1e2
, 1e3
, 1e4
, 1e5
, 1e6
, 1e7
, 1e8
, 1e9
,
1170 1e10
, 1e11
, 1e12
, 1e13
, 1e14
, 1e15
, 1e16
, 1e17
, 1e18
, 1e19
,
1178 static CONST
double bigtens
[] = { 1e16
, 1e32
, 1e64
, 1e128
, 1e256
};
1179 static CONST
double tinytens
[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 };
1183 static CONST
double bigtens
[] = { 1e16
, 1e32
, 1e64
};
1184 static CONST
double tinytens
[] = { 1e-16, 1e-32, 1e-64 };
1187 static CONST
double bigtens
[] = { 1e16
, 1e32
};
1188 static CONST
double tinytens
[] = { 1e-16, 1e-32 };
1194 static int quorem(Bigint
*b
, Bigint
*S
)
1199 ULong
*bx
, *bxe
, *sx
, *sxe
;
1212 q
= *bxe
/ (*sxe
+ 1); /* ensure q <= true quotient */
1219 ys
= (si
& 0xffff) * q
+ carry
;
1220 zs
= (si
>> 16) * q
+ (ys
>> 16);
1222 y
= (*bx
& 0xffff) - (ys
& 0xffff) + borrow
;
1224 Sign_Extend(borrow
, y
);
1225 z
= (*bx
>> 16) - (zs
& 0xffff) + borrow
;
1227 Sign_Extend(borrow
, z
);
1230 ys
= *sx
++ * q
+ carry
;
1232 y
= *bx
- (ys
& 0xffff) + borrow
;
1234 Sign_Extend(borrow
, y
);
1241 while(--bxe
> bx
&& !*bxe
)
1246 if (cmp(b
, S
) >= 0) {
1255 ys
= (si
& 0xffff) + carry
;
1256 zs
= (si
>> 16) + (ys
>> 16);
1258 y
= (*bx
& 0xffff) - (ys
& 0xffff) + borrow
;
1260 Sign_Extend(borrow
, y
);
1261 z
= (*bx
>> 16) - (zs
& 0xffff) + borrow
;
1263 Sign_Extend(borrow
, z
);
1268 y
= *bx
- (ys
& 0xffff) + borrow
;
1270 Sign_Extend(borrow
, y
);
1278 while(--bxe
> bx
&& !*bxe
)
1286 void destroy_freelist(Bigint
** freelist
)
1291 for (i
= 0; i
<= Kmax
; i
++) {
1292 Bigint
**listp
= &freelist
[i
];
1293 while ((tmp
= *listp
) != nullptr) {
1297 freelist
[i
] = nullptr;
1303 void zend_freedtoa(char *s
)
1305 Bigint
*b
= (Bigint
*)((int *)s
- 1);
1306 b
->maxwds
= 1 << (b
->k
= *(int*)b
);
1310 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
1312 * Inspired by "How to Print Floating-Point Numbers Accurately" by
1313 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
1316 * 1. Rather than iterating, we use a simple numeric overestimate
1317 * to determine k = floor(log10(d)). We scale relevant
1318 * quantities using O(log2(k)) rather than O(k) multiplications.
1319 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
1320 * try to generate digits strictly left to right. Instead, we
1321 * compute with fewer bits and propagate the carry if necessary
1322 * when rounding the final digit up. This is often faster.
1323 * 3. Under the assumption that input will be rounded nearest,
1324 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
1325 * That is, we allow equality in stopping tests when the
1326 * round-nearest rule will give the same floating-point value
1327 * as would satisfaction of the stopping test with strict
1329 * 4. We remove common factors of powers of 2 from relevant
1331 * 5. When converting floating-point integers less than 1e16,
1332 * we use floating-point arithmetic rather than resorting
1333 * to multiple-precision integers.
1334 * 6. When asked to produce fewer than 15 digits, we first try
1335 * to get by with floating-point arithmetic; we resort to
1336 * multiple-precision integer arithmetic only if we cannot
1337 * guarantee that the floating-point calculation has given
1338 * the correctly rounded result. For k requested digits and
1339 * "uniformly" distributed input, the probability is
1340 * something like 10^(k-15) that we must resort to the Long
1344 char * zend_dtoa(double _d
, int mode
, int ndigits
, int *decpt
, int *sign
, char **rve
)
1346 /* Arguments ndigits, decpt, sign are similar to those
1347 of ecvt and fcvt; trailing zeros are suppressed from
1348 the returned string. If not null, *rve is set to point
1349 to the end of the return value. If d is +-Infinity or NaN,
1350 then *decpt is set to 9999.
1353 0 ==> shortest string that yields d when read in
1354 and rounded to nearest.
1355 1 ==> like 0, but with Steele & White stopping rule;
1356 e.g. with IEEE P754 arithmetic , mode 0 gives
1357 1e23 whereas mode 1 gives 9.999999999999999e22.
1358 2 ==> max(1,ndigits) significant digits. This gives a
1359 return value similar to that of ecvt, except
1360 that trailing zeros are suppressed.
1361 3 ==> through ndigits past the decimal point. This
1362 gives a return value similar to that from fcvt,
1363 except that trailing zeros are suppressed, and
1364 ndigits can be negative.
1365 4-9 should give the same return values as 2-3, i.e.,
1366 4 <= mode <= 9 ==> same return as mode
1367 2 + (mode & 1). These modes are mainly for
1368 debugging; often they run slower but sometimes
1369 faster than modes 2-3.
1370 4,5,8,9 ==> left-to-right digit generation.
1371 6-9 ==> don't try fast floating-point estimate
1374 Values of mode other than 0-9 are treated as mode 0.
1376 Sufficient space is allocated to the return value
1377 to hold the suppressed trailing zeros.
1380 int bbits
, b2
, b5
, be
, dig
, i
, ieps
, ilim
= 0, ilim0
, ilim1
= 0,
1381 j
, j1
, k
, k0
, k_check
, leftright
, m2
, m5
, s2
, s5
,
1382 spec_case
= 0, try_quick
;
1384 #ifndef Sudden_Underflow
1388 Bigint
*b
, *b1
, *delta
, *mlo
= 0, *mhi
, *S
, *tmp
;
1395 if (word0(d
) & Sign_bit
) {
1396 /* set sign for everything, including 0's and NaNs */
1398 word0(d
) &= ~Sign_bit
; /* clear sign bit */
1403 #if defined(IEEE_Arith) + defined(VAX)
1405 if ((word0(d
) & Exp_mask
) == Exp_mask
)
1407 if (word0(d
) == 0x8000)
1410 /* Infinity or NaN */
1413 if (!word1(d
) && !(word0(d
) & 0xfffff))
1414 return nrv_alloc("Infinity", rve
, 8);
1416 return nrv_alloc("NaN", rve
, 3);
1420 value(d
) += 0; /* normalize */
1424 return nrv_alloc("0", rve
, 1);
1427 b
= d2b(value(d
), &be
, &bbits
);
1428 #ifdef Sudden_Underflow
1429 i
= (int)(word0(d
) >> Exp_shift1
& (Exp_mask
>>Exp_shift1
));
1431 if ((i
= (int)(word0(d
) >> Exp_shift1
& (Exp_mask
>>Exp_shift1
)))) {
1433 value(d2
) = value(d
);
1434 word0(d2
) &= Frac_mask1
;
1435 word0(d2
) |= Exp_11
;
1437 if (j
= 11 - hi0bits(word0(d2
) & Frac_mask
))
1438 value(d2
) /= 1 << j
;
1441 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
1442 * log10(x) = log(x) / log(10)
1443 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
1444 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
1446 * This suggests computing an approximation k to log10(d) by
1448 * k = (i - Bias)*0.301029995663981
1449 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
1451 * We want k to be too large rather than too small.
1452 * The error in the first-order Taylor series approximation
1453 * is in our favor, so we just round up the constant enough
1454 * to compensate for any error in the multiplication of
1455 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
1456 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
1457 * adding 1e-13 to the constant term more than suffices.
1458 * Hence we adjust the constant term to 0.1760912590558.
1459 * (We could get a more accurate k by invoking log10,
1460 * but this is probably not worthwhile.)
1468 #ifndef Sudden_Underflow
1472 /* d is denormalized */
1474 i
= bbits
+ be
+ (Bias
+ (P
-1) - 1);
1475 x
= i
> 32 ? (word0(d
) << (64 - i
)) | (word1(d
) >> (i
- 32))
1476 : (word1(d
) << (32 - i
));
1478 word0(d2
) -= 31*Exp_msk1
; /* adjust exponent */
1479 i
-= (Bias
+ (P
-1) - 1) + 1;
1483 ds
= (value(d2
)-1.5)*0.289529654602168 + 0.1760912590558 + i
*0.301029995663981;
1485 if (ds
< 0. && ds
!= k
)
1486 k
--; /* want k = floor(ds) */
1488 if (k
>= 0 && k
<= Ten_pmax
) {
1489 if (value(d
) < tens
[k
])
1512 if (mode
< 0 || mode
> 9)
1533 ilim
= ilim1
= i
= ndigits
;
1539 i
= ndigits
+ k
+ 1;
1545 s
= s0
= rv_alloc(i
);
1547 if (ilim
>= 0 && ilim
<= Quick_max
&& try_quick
) {
1549 /* Try to get by with floating-point arithmetic. */
1552 value(d2
) = value(d
);
1555 ieps
= 2; /* conservative */
1560 /* prevent overflows */
1562 value(d
) /= bigtens
[n_bigtens
-1];
1565 for(; j
; j
>>= 1, i
++)
1572 else if ((j1
= -k
)) {
1573 value(d
) *= tens
[j1
& 0xf];
1574 for(j
= j1
>> 4; j
; j
>>= 1, i
++)
1577 value(d
) *= bigtens
[i
];
1580 if (k_check
&& value(d
) < 1. && ilim
> 0) {
1588 value(eps
) = ieps
*value(d
) + 7.;
1589 word0(eps
) -= (P
-1)*Exp_msk1
;
1593 if (value(d
) > value(eps
))
1595 if (value(d
) < -value(eps
))
1599 #ifndef No_leftright
1601 /* Use Steele & White method of only
1602 * generating digits needed.
1604 value(eps
) = 0.5/tens
[ilim
-1] - value(eps
);
1606 L
= (int32_t)value(d
);
1608 *s
++ = '0' + (int)L
;
1609 if (value(d
) < value(eps
))
1611 if (1. - value(d
) < value(eps
))
1621 /* Generate ilim digits, then fix them up. */
1622 value(eps
) *= tens
[ilim
-1];
1623 for(i
= 1;; i
++, value(d
) *= 10.) {
1624 L
= (int32_t)value(d
);
1626 *s
++ = '0' + (int)L
;
1628 if (value(d
) > 0.5 + value(eps
))
1630 else if (value(d
) < 0.5 - value(eps
)) {
1631 /* cut ALL traling zeros only if the number of chars is greater than precision
1632 * otherwise cut only extra zeros
1635 while(*--s
== '0' && (s
- s0
) > k
);
1645 #ifndef No_leftright
1650 value(d
) = value(d2
);
1655 /* Do we have a "small" integer? */
1657 if (be
>= 0 && k
<= Int_max
) {
1660 if (ndigits
< 0 && ilim
<= 0) {
1662 if (ilim
< 0 || value(d
) <= 5*ds
)
1667 L
= (int32_t)(value(d
) / ds
);
1669 #ifdef Check_FLT_ROUNDS
1670 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
1676 *s
++ = '0' + (int)L
;
1678 value(d
) += value(d
);
1679 if (value(d
) > ds
|| (value(d
) == ds
&& (L
& 1))) {
1691 if (!(value(d
) *= 10.))
1703 #ifndef Sudden_Underflow
1704 denorm
? be
+ (Bias
+ (P
-1) - 1 + 1) :
1707 1 + 4*P
- 3 - bbits
+ ((bbits
+ be
- 1) & 3);
1721 if ((i
= ilim
) < 0) {
1730 if (m2
> 0 && s2
> 0) {
1731 i
= m2
< s2
? m2
: s2
;
1739 mhi
= pow5mult(mhi
, m5
);
1744 if ((j
= b5
- m5
)) {
1748 b
= pow5mult(b
, b5
);
1753 S
= pow5mult(S
, s5
);
1754 /* Check for special case that d is a normalized power of 2. */
1757 if (!word1(d
) && !(word0(d
) & Bndry_mask
)
1758 #ifndef Sudden_Underflow
1759 && word0(d
) & Exp_mask
1762 /* The special case */
1771 /* Arrange for convenient computation of quotients:
1772 * shift left if necessary so divisor has 4 leading 0 bits.
1774 * Perhaps we should just compute leading 28 bits of S once
1775 * and for all and pass them and a shift to quorem, so it
1776 * can do shifts and ors to compute the numerator for q.
1779 if ((i
= ((s5
? 32 - hi0bits(S
->x
[S
->wds
-1]) : 1) + s2
) & 0x1f))
1782 if ((i
= ((s5
? 32 - hi0bits(S
->x
[S
->wds
-1]) : 1) + s2
) & 0xf))
1804 b
= multadd(b
, 10, 0); /* we botched the k estimate */
1806 mhi
= multadd(mhi
, 10, 0);
1810 if (ilim
<= 0 && mode
> 2) {
1811 if (ilim
< 0 || cmp(b
,S
= multadd(S
,5,0)) <= 0) {
1812 /* no digits, fcvt style */
1824 mhi
= lshift(mhi
, m2
);
1826 /* Compute mlo -- check for special case
1827 * that d is a normalized power of 2.
1832 mhi
= Balloc(mhi
->k
);
1834 mhi
= lshift(mhi
, Log2P
);
1838 dig
= quorem(b
,S
) + '0';
1839 /* Do we yet have the shortest decimal string
1840 * that will round to d?
1843 delta
= diff(S
, mhi
);
1844 j1
= delta
->sign
? 1 : cmp(b
, delta
);
1846 #ifndef ROUND_BIASED
1847 if (j1
== 0 && !mode
&& !(word1(d
) & 1)) {
1856 if (j
< 0 || (j
== 0 && !mode
1857 #ifndef ROUND_BIASED
1864 if ((j1
> 0 || (j1
== 0 && (dig
& 1)))
1872 if (dig
== '9') { /* possible if i == 1 */
1883 b
= multadd(b
, 10, 0);
1885 mlo
= mhi
= multadd(mhi
, 10, 0);
1887 mlo
= multadd(mlo
, 10, 0);
1888 mhi
= multadd(mhi
, 10, 0);
1894 *s
++ = dig
= quorem(b
,S
) + '0';
1897 b
= multadd(b
, 10, 0);
1900 /* Round off last digit */
1904 if (j
> 0 || (j
== 0 && (dig
& 1))) {
1921 if (mlo
&& mlo
!= mhi
)
1928 Bigint
*&p5s
= s_bigint_data
->p5s
;
1938 if (s
== s0
) { /* don't return empty string */
1949 double zend_strtod (CONST
char *s00
, const char **se
)
1951 int bb2
, bb5
, bbe
, bd2
, bd5
, bbbits
, bs2
, c
, dsign
,
1952 e
, e1
, esign
, i
, j
, k
, nd
, nd0
, nf
, nz
, nz0
, sign
;
1953 CONST
char *s
, *s0
, *s1
;
1954 double aadj
, aadj1
, adj
;
1958 Bigint
*bb
= 0, *bb1
, *bd
= 0, *bd0
, *bs
= 0, *delta
= 0, *tmp
;
1961 CONST
char decimal_point
= '.';
1963 sign
= nz0
= nz
= 0;
1967 for(s
= s00
; isspace((unsigned char) *s
); s
++)
1973 } else if (*s
== '+') {
1984 while(*++s
== '0') ;
1990 for(nd
= nf
= 0; (c
= *s
) >= '0' && c
<= '9'; nd
++, s
++)
1996 if (c
== decimal_point
) {
1999 for(; c
== '0'; c
= *++s
)
2001 if (c
> '0' && c
<= '9') {
2009 for(; c
>= '0' && c
<= '9'; c
= *++s
) {
2014 for(i
= 1; i
< nz
; i
++)
2017 else if (nd
<= DBL_DIG
+ 1)
2021 else if (nd
<= DBL_DIG
+ 1)
2029 if (c
== 'e' || c
== 'E') {
2030 if (!nd
&& !nz
&& !nz0
) {
2042 if (c
>= '0' && c
<= '9') {
2045 if (c
> '0' && c
<= '9') {
2048 while((c
= *++s
) >= '0' && c
<= '9')
2050 if (s
- s1
> 8 || L
> 19999)
2051 /* Avoid confusion from exponents
2052 * so large that e might overflow.
2054 e
= 19999; /* safe for 16 bit ints */
2073 /* Now we have nd0 digits, starting at s0, followed by a
2074 * decimal point, followed by nd-nd0 digits. The number we're
2075 * after is the integer represented by those digits times
2080 k
= nd
< DBL_DIG
+ 1 ? nd
: DBL_DIG
+ 1;
2083 value(rv
) = tens
[k
- 9] * value(rv
) + z
;
2086 #ifndef RND_PRODQUOT
2093 if (e
<= Ten_pmax
) {
2095 goto vax_ovfl_check
;
2097 /* value(rv) = */ rounded_product(value(rv
),
2103 if (e
<= Ten_pmax
+ i
) {
2104 /* A fancier test would sometimes let us do
2105 * this for larger i values.
2108 value(rv
) *= tens
[i
];
2110 /* VAX exponent range is so narrow we must
2111 * worry about overflow here...
2114 word0(rv
) -= P
*Exp_msk1
;
2115 /* value(rv) = */ rounded_product(value(rv
),
2117 if ((word0(rv
) & Exp_mask
)
2118 > Exp_msk1
*(DBL_MAX_EXP
+Bias
-1-P
))
2120 word0(rv
) += P
*Exp_msk1
;
2122 /* value(rv) = */ rounded_product(value(rv
),
2128 #ifndef Inaccurate_Divide
2129 else if (e
>= -Ten_pmax
) {
2130 /* value(rv) = */ rounded_quotient(value(rv
),
2138 /* Get starting approximation = rv * 10**e1 */
2142 value(rv
) *= tens
[i
];
2144 if (e1
> DBL_MAX_10_EXP
) {
2148 value(rv
) = HUGE_VAL
;
2150 /* Can't trust HUGE_VAL */
2152 word0(rv
) = Exp_mask
;
2164 for(j
= 0; e1
> 1; j
++, e1
>>= 1)
2166 value(rv
) *= bigtens
[j
];
2167 /* The last multiplication could overflow. */
2168 word0(rv
) -= P
*Exp_msk1
;
2169 value(rv
) *= bigtens
[j
];
2170 if ((z
= word0(rv
) & Exp_mask
)
2171 > Exp_msk1
*(DBL_MAX_EXP
+Bias
-P
))
2173 if (z
> Exp_msk1
*(DBL_MAX_EXP
+Bias
-1-P
)) {
2174 /* set to largest number */
2175 /* (Can't trust DBL_MAX) */
2180 word0(rv
) += P
*Exp_msk1
;
2188 value(rv
) /= tens
[i
];
2191 if (e1
>= 1 << n_bigtens
)
2193 for(j
= 0; e1
> 1; j
++, e1
>>= 1)
2195 value(rv
) *= tinytens
[j
];
2196 /* The last multiplication could underflow. */
2197 value(rv0
) = value(rv
);
2198 value(rv
) *= tinytens
[j
];
2200 value(rv
) = 2.*value(rv0
);
2201 value(rv
) *= tinytens
[j
];
2212 /* The refinement below will clean
2213 * this approximation up.
2219 /* Now the hard part -- adjusting rv to the correct value.*/
2221 /* Put digits into bd: true value = bd * 10^e */
2223 bd0
= s2b(s0
, nd0
, nd
, y
);
2226 bd
= Balloc(bd0
->k
);
2228 bb
= d2b(value(rv
), &bbe
, &bbbits
); /* rv = bb * 2^bbe */
2244 #ifdef Sudden_Underflow
2246 j
= 1 + 4*P
- 3 - bbbits
+ ((bbe
+ bbbits
- 1) & 3);
2251 i
= bbe
+ bbbits
- 1; /* logb(rv) */
2252 if (i
< Emin
) /* denormal */
2259 i
= bb2
< bd2
? bb2
: bd2
;
2268 bs
= pow5mult(bs
, bb5
);
2274 bb
= lshift(bb
, bb2
);
2276 bd
= pow5mult(bd
, bd5
);
2278 bd
= lshift(bd
, bd2
);
2280 bs
= lshift(bs
, bs2
);
2281 delta
= diff(bb
, bd
);
2282 dsign
= delta
->sign
;
2286 /* Error is less than half an ulp -- check for
2287 * special case of mantissa a power of two.
2289 if (dsign
|| word1(rv
) || word0(rv
) & Bndry_mask
)
2291 delta
= lshift(delta
,Log2P
);
2292 if (cmp(delta
, bs
) > 0)
2297 /* exactly half-way between */
2299 if ((word0(rv
) & Bndry_mask1
) == Bndry_mask1
2300 && word1(rv
) == 0xffffffff) {
2301 /*boundary case -- increment exponent*/
2302 word0(rv
) = (word0(rv
) & Exp_mask
)
2312 else if (!(word0(rv
) & Bndry_mask
) && !word1(rv
)) {
2314 /* boundary case -- decrement exponent */
2315 #ifdef Sudden_Underflow
2316 L
= word0(rv
) & Exp_mask
;
2325 L
= (word0(rv
) & Exp_mask
) - Exp_msk1
;
2327 word0(rv
) = L
| Bndry_mask1
;
2328 word1(rv
) = 0xffffffff;
2335 #ifndef ROUND_BIASED
2336 if (!(word1(rv
) & LSB
))
2340 value(rv
) += ulp(value(rv
));
2341 #ifndef ROUND_BIASED
2343 value(rv
) -= ulp(value(rv
));
2344 #ifndef Sudden_Underflow
2352 if ((aadj
= ratio(delta
, bs
)) <= 2.) {
2355 else if (word1(rv
) || word0(rv
) & Bndry_mask
) {
2356 #ifndef Sudden_Underflow
2357 if (word1(rv
) == Tiny1
&& !word0(rv
))
2364 /* special case -- power of FLT_RADIX to be */
2365 /* rounded down... */
2367 if (aadj
< 2./FLT_RADIX
)
2368 aadj
= 1./FLT_RADIX
;
2376 aadj1
= dsign
? aadj
: -aadj
;
2377 #ifdef Check_FLT_ROUNDS
2378 switch(FLT_ROUNDS
) {
2379 case 2: /* towards +infinity */
2382 case 0: /* towards 0 */
2383 case 3: /* towards -infinity */
2387 if (FLT_ROUNDS
== 0)
2391 y
= word0(rv
) & Exp_mask
;
2393 /* Check for overflow */
2395 if (y
== Exp_msk1
*(DBL_MAX_EXP
+Bias
-1)) {
2396 value(rv0
) = value(rv
);
2397 word0(rv
) -= P
*Exp_msk1
;
2398 adj
= aadj1
* ulp(value(rv
));
2400 if ((word0(rv
) & Exp_mask
) >=
2401 Exp_msk1
*(DBL_MAX_EXP
+Bias
-P
)) {
2402 if (word0(rv0
) == Big0
&& word1(rv0
) == Big1
)
2409 word0(rv
) += P
*Exp_msk1
;
2412 #ifdef Sudden_Underflow
2413 if ((word0(rv
) & Exp_mask
) <= P
*Exp_msk1
) {
2414 value(rv0
) = value(rv
);
2415 word0(rv
) += P
*Exp_msk1
;
2416 adj
= aadj1
* ulp(value(rv
));
2419 if ((word0(rv
) & Exp_mask
) < P
*Exp_msk1
)
2421 if ((word0(rv
) & Exp_mask
) <= P
*Exp_msk1
)
2424 if (word0(rv0
) == Tiny0
2425 && word1(rv0
) == Tiny1
)
2432 word0(rv
) -= P
*Exp_msk1
;
2435 adj
= aadj1
* ulp(value(rv
));
2439 /* Compute adj so that the IEEE rounding rules will
2440 * correctly round rv + adj in some half-way cases.
2441 * If rv * ulp(rv) is denormalized (i.e.,
2442 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
2443 * trouble from bits lost to denormalization;
2444 * example: 1.2e-307 .
2446 if (y
<= (P
-1)*Exp_msk1
&& aadj
>= 1.) {
2447 aadj1
= (double)(int)(aadj
+ 0.5);
2451 adj
= aadj1
* ulp(value(rv
));
2455 z
= word0(rv
) & Exp_mask
;
2457 /* Can we stop now? */
2460 /* The tolerances below are conservative. */
2461 if (dsign
|| word1(rv
) || word0(rv
) & Bndry_mask
) {
2462 if (aadj
< .4999999 || aadj
> .5000001)
2465 else if (aadj
< .4999999/FLT_RADIX
)
2483 result
= sign
? -value(rv
) : value(rv
);
2485 if (s_bigint_data
.isNull()) {
2489 Bigint
*&p5s
= s_bigint_data
->p5s
;
2499 double zend_hex_strtod(const char *str
, const char **endptr
)
2501 const char *s
= str
;
2506 if (*s
== '0' && (s
[1] == 'x' || s
[1] == 'X')) {
2510 while ((c
= *s
++)) {
2511 if (c
>= '0' && c
<= '9') {
2513 } else if (c
>= 'A' && c
<= 'F') {
2515 } else if (c
>= 'a' && c
<= 'f') {
2522 value
= value
* 16 + c
;
2525 if (endptr
!= nullptr) {
2526 *endptr
= (char *)(any
? s
- 1 : str
);
2532 double zend_oct_strtod(const char *str
, const char **endptr
)
2534 const char *s
= str
;
2539 /* skip leading zero */
2542 while ((c
= *s
++)) {
2544 /* break and return the current value if the number is not well-formed
2545 * that's what Linux strtol() does
2549 // Note parentheses: convert digit into integer before adding as a double
2550 // in order to avoid accumulating floating point inaccuracies
2551 value
= value
* 8 + (c
- '0');
2555 if (endptr
!= nullptr) {
2556 *endptr
= (char *)(any
? s
- 1 : str
);
2562 double zend_bin_strtod(const char *str
, const char **endptr
)
2564 const char *s
= str
;
2569 if (strlen(str
) < 2) {
2574 if ('0' == *s
&& ('b' == s
[1] || 'B' == s
[1])) {
2578 while ((c
= *s
++)) {
2580 * Verify the validity of the current character as a base-2 digit. In
2581 * the event that an invalid digit is found, halt the conversion and
2582 * return the portion which has been converted thus far.
2584 if ('0' == c
|| '1' == c
)
2585 // Note parentheses: convert digit into integer before adding as a double
2586 // in order to avoid accumulating floating point inaccuracies
2587 value
= value
* 2 + (c
- '0');
2595 * As with many strtoX implementations, should the subject sequence be
2596 * empty or not well-formed, no conversion is performed and the original
2597 * value of str is stored in *endptr, provided that endptr is not a null
2600 if (nullptr != endptr
) {
2601 *endptr
= (char *)(any
? s
- 1 : str
);
2607 ///////////////////////////////////////////////////////////////////////////////