From 1ca7c026d1cd8ea06c528ffd2cf81cc7dbbc9627 Mon Sep 17 00:00:00 2001 From: Paul Tarjan Date: Fri, 14 Jun 2013 18:10:19 -0700 Subject: [PATCH] Revert double thing Some tests are failing now that didn't fail in contbuild. Revert until I can look at it. --- hphp/runtime/base/string_data.cpp | 2 +- hphp/runtime/base/variable_unserializer.cpp | 2 +- hphp/runtime/base/zend/zend_collator.cpp | 2 +- hphp/runtime/base/zend/zend_functions.cpp | 2 +- hphp/runtime/base/zend/zend_strtod.cpp | 719 ++++++++++++++++++++ hphp/runtime/base/zend/zend_strtod.h | 1 + 6 files changed, 724 insertions(+), 4 deletions(-) diff --git a/hphp/runtime/base/string_data.cpp b/hphp/runtime/base/string_data.cpp index 28c361531..59ebe8eb5 100644 --- a/hphp/runtime/base/string_data.cpp +++ b/hphp/runtime/base/string_data.cpp @@ -824,7 +824,7 @@ int64_t StringData::toInt64(int base /* = 10 */) const { double StringData::toDouble() const { StringSlice s = slice(); - if (s.len) return strtod(s.ptr, nullptr); + if (s.len) return zend_strtod(s.ptr, nullptr); return 0; } diff --git a/hphp/runtime/base/variable_unserializer.cpp b/hphp/runtime/base/variable_unserializer.cpp index 02aed7eab..3cbd1300c 100644 --- a/hphp/runtime/base/variable_unserializer.cpp +++ b/hphp/runtime/base/variable_unserializer.cpp @@ -46,7 +46,7 @@ int64_t VariableUnserializer::readInt() { double VariableUnserializer::readDouble() { check(); char *newBuf; - double r = strtod(m_buf, &newBuf); + double r = zend_strtod(m_buf, &newBuf); m_buf = newBuf; return r; } diff --git a/hphp/runtime/base/zend/zend_collator.cpp b/hphp/runtime/base/zend/zend_collator.cpp index 5c47c2118..6abaccf66 100644 --- a/hphp/runtime/base/zend/zend_collator.cpp +++ b/hphp/runtime/base/zend/zend_collator.cpp @@ -98,7 +98,7 @@ static double collator_u_strtod(const UChar *nptr, UChar **endptr) { } *bufpos = '\0'; - value = strtod(numbuf, nullptr); + value = zend_strtod(numbuf, nullptr); if (numbuf != buf) { free(numbuf); diff --git a/hphp/runtime/base/zend/zend_functions.cpp b/hphp/runtime/base/zend/zend_functions.cpp index bceffdf22..711da4706 100644 --- a/hphp/runtime/base/zend/zend_functions.cpp +++ b/hphp/runtime/base/zend/zend_functions.cpp @@ -149,7 +149,7 @@ DataType is_numeric_string(const char *str, int length, int64_t *lval, } else if (!(digits < SIZEOF_LONG * 2 || (digits == SIZEOF_LONG * 2 && ptr[-digits] <= '7'))) { if (dval) { - local_dval = strtod(str, (char **)&ptr); + local_dval = zend_strtod(str, (char **)&ptr); } type = KindOfDouble; } diff --git a/hphp/runtime/base/zend/zend_strtod.cpp b/hphp/runtime/base/zend/zend_strtod.cpp index 94748078a..922bb1d99 100644 --- a/hphp/runtime/base/zend/zend_strtod.cpp +++ b/hphp/runtime/base/zend/zend_strtod.cpp @@ -257,6 +257,7 @@ extern void *MALLOC(size_t); #define Emin (-1022) #define Exp_1 0x3ff00000 #define Exp_11 0x3ff00000 +#define Ebits 11 #define Frac_mask 0xfffff #define Frac_mask1 0xfffff #define Ten_pmax 22 @@ -284,6 +285,7 @@ extern void *MALLOC(size_t); #define Bias 65 #define Exp_1 0x41000000 #define Exp_11 0x41000000 +#define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */ #define Frac_mask 0xffffff #define Frac_mask1 0xffffff #define Bletch 4 @@ -307,6 +309,7 @@ extern void *MALLOC(size_t); #define Bias 129 #define Exp_1 0x40800000 #define Exp_11 0x4080 +#define Ebits 8 #define Frac_mask 0x7fffff #define Frac_mask1 0xffff007f #define Ten_pmax 24 @@ -660,6 +663,39 @@ static Bigint * mult(Bigint *a, Bigint *b) return c; } +static Bigint * s2b (CONST char *s, int nd0, int nd, ULong y9) +{ + Bigint *b; + int i, k; + Long x, y; + + x = (nd + 8) / 9; + for(k = 0, y = 1; x > y; y <<= 1, k++) ; +#ifdef Pack_32 + b = Balloc(k); + b->x[0] = y9; + b->wds = 1; +#else + b = Balloc(k+1); + b->x[0] = y9 & 0xffff; + b->wds = (b->x[1] = y9 >> 16) ? 2 : 1; +#endif + + i = 9; + if (9 < nd0) { + s += 9; + do b = multadd(b, 10, *s++ - '0'); + while(++i < nd0); + s++; + } else { + s += 10; + } + for(; i < nd; i++) { + b = multadd(b, 10, *s++ - '0'); + } + return b; +} + static Bigint * pow5mult(Bigint *b, int k) { Bigint *b1, *p5, *p51; @@ -856,6 +892,108 @@ static Bigint * diff(Bigint *a, Bigint *b) return c; } +static double ulp (double _x) +{ + _double x; + register Long L; + _double a; + + value(x) = _x; + L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1; +#ifndef Sudden_Underflow + if (L > 0) { +#endif +#ifdef IBM + L |= Exp_msk1 >> 4; +#endif + word0(a) = L; + word1(a) = 0; +#ifndef Sudden_Underflow + } + else { + L = -L >> Exp_shift; + if (L < Exp_shift) { + word0(a) = 0x80000 >> L; + word1(a) = 0; + } + else { + word0(a) = 0; + L -= Exp_shift; + word1(a) = L >= 31 ? 1 : 1 << (31 - L); + } + } +#endif + return value(a); +} + +static double +b2d +#ifdef KR_headers +(a, e) Bigint *a; int *e; +#else +(Bigint *a, int *e) +#endif +{ + ULong *xa, *xa0, w, y, z; + int k; + _double d; +#ifdef VAX + ULong d0, d1; +#else +#define d0 word0(d) +#define d1 word1(d) +#endif + + xa0 = a->x; + xa = xa0 + a->wds; + y = *--xa; + k = hi0bits(y); + *e = 32 - k; +#ifdef Pack_32 + if (k < Ebits) { + d0 = Exp_1 | y >> (Ebits - k); + w = xa > xa0 ? *--xa : 0; + d1 = y << ((32-Ebits) + k) | w >> (Ebits - k); + goto ret_d; + } + z = xa > xa0 ? *--xa : 0; + if (k -= Ebits) { + d0 = Exp_1 | y << k | z >> (32 - k); + y = xa > xa0 ? *--xa : 0; + d1 = z << k | y >> (32 - k); + } + else { + d0 = Exp_1 | y; + d1 = z; + } +#else + if (k < Ebits + 16) { + z = xa > xa0 ? *--xa : 0; + d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k; + w = xa > xa0 ? *--xa : 0; + y = xa > xa0 ? *--xa : 0; + d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k; + goto ret_d; + } + z = xa > xa0 ? *--xa : 0; + w = xa > xa0 ? *--xa : 0; + k -= Ebits + 16; + d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k; + y = xa > xa0 ? *--xa : 0; + d1 = w << k + 16 | y << k; +#endif +ret_d: +#ifdef VAX + word0(d) = d0 >> 16 | d0 << 16; + word1(d) = d1 >> 16 | d1 << 16; +#else +#undef d0 +#undef d1 +#endif + return value(d); +} + + static Bigint * d2b(double _d, int *e, int *bits) { Bigint *b; @@ -972,6 +1110,41 @@ static Bigint * d2b(double _d, int *e, int *bits) #undef d1 +static double ratio (Bigint *a, Bigint *b) +{ + _double da, db; + int k, ka, kb; + + value(da) = b2d(a, &ka); + value(db) = b2d(b, &kb); +#ifdef Pack_32 + k = ka - kb + 32*(a->wds - b->wds); +#else + k = ka - kb + 16*(a->wds - b->wds); +#endif +#ifdef IBM + if (k > 0) { + word0(da) += (k >> 2)*Exp_msk1; + if (k &= 3) { + da *= 1 << k; + } + } else { + k = -k; + word0(db) += (k >> 2)*Exp_msk1; + if (k &= 3) + db *= 1 << k; + } +#else + if (k > 0) { + word0(da) += k*Exp_msk1; + } else { + k = -k; + word0(db) += k*Exp_msk1; + } +#endif + return value(da) / value(db); +} + static CONST double tens[] = { 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, @@ -1755,6 +1928,552 @@ ret1: return s0; } +double zend_strtod (CONST char *s00, char **se) +{ + int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign, + e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign; + CONST char *s, *s0, *s1; + double aadj, aadj1, adj; + _double rv, rv0; + Long L; + ULong y, z; + Bigint *bb = 0, *bb1, *bd = 0, *bd0, *bs = 0, *delta = 0, *tmp; + double result; + + CONST char decimal_point = '.'; + + sign = nz0 = nz = 0; + value(rv) = 0.; + + + for(s = s00; isspace((unsigned char) *s); s++) + ; + + if (*s == '-') { + sign = 1; + s++; + } else if (*s == '+') { + s++; + } + + if (*s == '\0') { + s = s00; + goto ret; + } + + if (*s == '0') { + nz0 = 1; + while(*++s == '0') ; + if (!*s) + goto ret; + } + s0 = s; + y = z = 0; + for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) + if (nd < 9) + y = 10*y + c - '0'; + else if (nd < 16) + z = 10*z + c - '0'; + nd0 = nd; + if (c == decimal_point) { + c = *++s; + if (!nd) { + for(; c == '0'; c = *++s) + nz++; + if (c > '0' && c <= '9') { + s0 = s; + nf += nz; + nz = 0; + goto have_dig; + } + goto dig_done; + } + for(; c >= '0' && c <= '9'; c = *++s) { +have_dig: + nz++; + if (c -= '0') { + nf += nz; + for(i = 1; i < nz; i++) + if (nd++ < 9) + y *= 10; + else if (nd <= DBL_DIG + 1) + z *= 10; + if (nd++ < 9) + y = 10*y + c; + else if (nd <= DBL_DIG + 1) + z = 10*z + c; + nz = 0; + } + } + } +dig_done: + e = 0; + if (c == 'e' || c == 'E') { + if (!nd && !nz && !nz0) { + s = s00; + goto ret; + } + s00 = s; + esign = 0; + switch(c = *++s) { + case '-': + esign = 1; + case '+': + c = *++s; + } + if (c >= '0' && c <= '9') { + while(c == '0') + c = *++s; + if (c > '0' && c <= '9') { + L = c - '0'; + s1 = s; + while((c = *++s) >= '0' && c <= '9') + L = 10*L + c - '0'; + if (s - s1 > 8 || L > 19999) + /* Avoid confusion from exponents + * so large that e might overflow. + */ + e = 19999; /* safe for 16 bit ints */ + else + e = (int)L; + if (esign) + e = -e; + } + else + e = 0; + } + else + s = s00; + } + if (!nd) { + if (!nz && !nz0) + s = s00; + goto ret; + } + e1 = e -= nf; + + /* Now we have nd0 digits, starting at s0, followed by a + * decimal point, followed by nd-nd0 digits. The number we're + * after is the integer represented by those digits times + * 10**e */ + + if (!nd0) + nd0 = nd; + k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; + value(rv) = y; + if (k > 9) + value(rv) = tens[k - 9] * value(rv) + z; + bd0 = 0; + if (nd <= DBL_DIG +#ifndef RND_PRODQUOT + && FLT_ROUNDS == 1 +#endif + ) { + if (!e) + goto ret; + if (e > 0) { + if (e <= Ten_pmax) { +#ifdef VAX + goto vax_ovfl_check; +#else + /* value(rv) = */ rounded_product(value(rv), + tens[e]); + goto ret; +#endif + } + i = DBL_DIG - nd; + if (e <= Ten_pmax + i) { + /* A fancier test would sometimes let us do + * this for larger i values. + */ + e -= i; + value(rv) *= tens[i]; +#ifdef VAX + /* VAX exponent range is so narrow we must + * worry about overflow here... + */ +vax_ovfl_check: + word0(rv) -= P*Exp_msk1; + /* value(rv) = */ rounded_product(value(rv), + tens[e]); + if ((word0(rv) & Exp_mask) + > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) + goto ovfl; + word0(rv) += P*Exp_msk1; +#else + /* value(rv) = */ rounded_product(value(rv), + tens[e]); +#endif + goto ret; + } + } +#ifndef Inaccurate_Divide + else if (e >= -Ten_pmax) { + /* value(rv) = */ rounded_quotient(value(rv), + tens[-e]); + goto ret; + } +#endif + } + e1 += nd - k; + + /* Get starting approximation = rv * 10**e1 */ + + if (e1 > 0) { + if ((i = e1 & 15)) + value(rv) *= tens[i]; + if (e1 &= ~15) { + if (e1 > DBL_MAX_10_EXP) { +ovfl: + errno = ERANGE; +#ifndef Bad_float_h + value(rv) = HUGE_VAL; +#else + /* Can't trust HUGE_VAL */ +#ifdef IEEE_Arith + word0(rv) = Exp_mask; + word1(rv) = 0; +#else + word0(rv) = Big0; + word1(rv) = Big1; +#endif +#endif + if (bd0) + goto retfree; + goto ret; + } + if (e1 >>= 4) { + for(j = 0; e1 > 1; j++, e1 >>= 1) + if (e1 & 1) + value(rv) *= bigtens[j]; + /* The last multiplication could overflow. */ + word0(rv) -= P*Exp_msk1; + value(rv) *= bigtens[j]; + if ((z = word0(rv) & Exp_mask) + > Exp_msk1*(DBL_MAX_EXP+Bias-P)) + goto ovfl; + if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) { + /* set to largest number */ + /* (Can't trust DBL_MAX) */ + word0(rv) = Big0; + word1(rv) = Big1; + } + else + word0(rv) += P*Exp_msk1; + } + + } + } + else if (e1 < 0) { + e1 = -e1; + if ((i = e1 & 15)) + value(rv) /= tens[i]; + if (e1 &= ~15) { + e1 >>= 4; + if (e1 >= 1 << n_bigtens) + goto undfl; + for(j = 0; e1 > 1; j++, e1 >>= 1) + if (e1 & 1) + value(rv) *= tinytens[j]; + /* The last multiplication could underflow. */ + value(rv0) = value(rv); + value(rv) *= tinytens[j]; + if (!value(rv)) { + value(rv) = 2.*value(rv0); + value(rv) *= tinytens[j]; + if (!value(rv)) { +undfl: + value(rv) = 0.; + errno = ERANGE; + if (bd0) + goto retfree; + goto ret; + } + word0(rv) = Tiny0; + word1(rv) = Tiny1; + /* The refinement below will clean + * this approximation up. + */ + } + } + } + + /* Now the hard part -- adjusting rv to the correct value.*/ + + /* Put digits into bd: true value = bd * 10^e */ + + bd0 = s2b(s0, nd0, nd, y); + + for(;;) { + bd = Balloc(bd0->k); + Bcopy(bd, bd0); + bb = d2b(value(rv), &bbe, &bbbits); /* rv = bb * 2^bbe */ + bs = i2b(1); + + if (e >= 0) { + bb2 = bb5 = 0; + bd2 = bd5 = e; + } + else { + bb2 = bb5 = -e; + bd2 = bd5 = 0; + } + if (bbe >= 0) + bb2 += bbe; + else + bd2 -= bbe; + bs2 = bb2; +#ifdef Sudden_Underflow +#ifdef IBM + j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3); +#else + j = P + 1 - bbbits; +#endif +#else + i = bbe + bbbits - 1; /* logb(rv) */ + if (i < Emin) /* denormal */ + j = bbe + (P-Emin); + else + j = P + 1 - bbbits; +#endif + bb2 += j; + bd2 += j; + i = bb2 < bd2 ? bb2 : bd2; + if (i > bs2) + i = bs2; + if (i > 0) { + bb2 -= i; + bd2 -= i; + bs2 -= i; + } + if (bb5 > 0) { + bs = pow5mult(bs, bb5); + bb1 = mult(bs, bb); + Bfree(bb); + bb = bb1; + } + if (bb2 > 0) + bb = lshift(bb, bb2); + if (bd5 > 0) + bd = pow5mult(bd, bd5); + if (bd2 > 0) + bd = lshift(bd, bd2); + if (bs2 > 0) + bs = lshift(bs, bs2); + delta = diff(bb, bd); + dsign = delta->sign; + delta->sign = 0; + i = cmp(delta, bs); + if (i < 0) { + /* Error is less than half an ulp -- check for + * special case of mantissa a power of two. + */ + if (dsign || word1(rv) || word0(rv) & Bndry_mask) + break; + delta = lshift(delta,Log2P); + if (cmp(delta, bs) > 0) + goto drop_down; + break; + } + if (i == 0) { + /* exactly half-way between */ + if (dsign) { + if ((word0(rv) & Bndry_mask1) == Bndry_mask1 + && word1(rv) == 0xffffffff) { + /*boundary case -- increment exponent*/ + word0(rv) = (word0(rv) & Exp_mask) + + Exp_msk1 +#ifdef IBM + | Exp_msk1 >> 4 +#endif + ; + word1(rv) = 0; + break; + } + } + else if (!(word0(rv) & Bndry_mask) && !word1(rv)) { +drop_down: + /* boundary case -- decrement exponent */ +#ifdef Sudden_Underflow + L = word0(rv) & Exp_mask; +#ifdef IBM + if (L < Exp_msk1) +#else + if (L <= Exp_msk1) +#endif + goto undfl; + L -= Exp_msk1; +#else + L = (word0(rv) & Exp_mask) - Exp_msk1; +#endif + word0(rv) = L | Bndry_mask1; + word1(rv) = 0xffffffff; +#ifdef IBM + goto cont; +#else + break; +#endif + } +#ifndef ROUND_BIASED + if (!(word1(rv) & LSB)) + break; +#endif + if (dsign) + value(rv) += ulp(value(rv)); +#ifndef ROUND_BIASED + else { + value(rv) -= ulp(value(rv)); +#ifndef Sudden_Underflow + if (!value(rv)) + goto undfl; +#endif + } +#endif + break; + } + if ((aadj = ratio(delta, bs)) <= 2.) { + if (dsign) + aadj = aadj1 = 1.; + else if (word1(rv) || word0(rv) & Bndry_mask) { +#ifndef Sudden_Underflow + if (word1(rv) == Tiny1 && !word0(rv)) + goto undfl; +#endif + aadj = 1.; + aadj1 = -1.; + } + else { + /* special case -- power of FLT_RADIX to be */ + /* rounded down... */ + + if (aadj < 2./FLT_RADIX) + aadj = 1./FLT_RADIX; + else + aadj *= 0.5; + aadj1 = -aadj; + } + } + else { + aadj *= 0.5; + aadj1 = dsign ? aadj : -aadj; +#ifdef Check_FLT_ROUNDS + switch(FLT_ROUNDS) { + case 2: /* towards +infinity */ + aadj1 -= 0.5; + break; + case 0: /* towards 0 */ + case 3: /* towards -infinity */ + aadj1 += 0.5; + } +#else + if (FLT_ROUNDS == 0) + aadj1 += 0.5; +#endif + } + y = word0(rv) & Exp_mask; + + /* Check for overflow */ + + if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) { + value(rv0) = value(rv); + word0(rv) -= P*Exp_msk1; + adj = aadj1 * ulp(value(rv)); + value(rv) += adj; + if ((word0(rv) & Exp_mask) >= + Exp_msk1*(DBL_MAX_EXP+Bias-P)) { + if (word0(rv0) == Big0 && word1(rv0) == Big1) + goto ovfl; + word0(rv) = Big0; + word1(rv) = Big1; + goto cont; + } + else + word0(rv) += P*Exp_msk1; + } + else { +#ifdef Sudden_Underflow + if ((word0(rv) & Exp_mask) <= P*Exp_msk1) { + value(rv0) = value(rv); + word0(rv) += P*Exp_msk1; + adj = aadj1 * ulp(value(rv)); + value(rv) += adj; +#ifdef IBM + if ((word0(rv) & Exp_mask) < P*Exp_msk1) +#else + if ((word0(rv) & Exp_mask) <= P*Exp_msk1) +#endif + { + if (word0(rv0) == Tiny0 + && word1(rv0) == Tiny1) + goto undfl; + word0(rv) = Tiny0; + word1(rv) = Tiny1; + goto cont; + } + else + word0(rv) -= P*Exp_msk1; + } + else { + adj = aadj1 * ulp(value(rv)); + value(rv) += adj; + } +#else + /* Compute adj so that the IEEE rounding rules will + * correctly round rv + adj in some half-way cases. + * If rv * ulp(rv) is denormalized (i.e., + * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid + * trouble from bits lost to denormalization; + * example: 1.2e-307 . + */ + if (y <= (P-1)*Exp_msk1 && aadj >= 1.) { + aadj1 = (double)(int)(aadj + 0.5); + if (!dsign) + aadj1 = -aadj1; + } + adj = aadj1 * ulp(value(rv)); + value(rv) += adj; +#endif + } + z = word0(rv) & Exp_mask; + if (y == z) { + /* Can we stop now? */ + L = (int32_t)aadj; + aadj -= L; + /* The tolerances below are conservative. */ + if (dsign || word1(rv) || word0(rv) & Bndry_mask) { + if (aadj < .4999999 || aadj > .5000001) + break; + } + else if (aadj < .4999999/FLT_RADIX) + break; + } +cont: + Bfree(bb); + Bfree(bd); + Bfree(bs); + Bfree(delta); + } +retfree: + Bfree(bb); + Bfree(bd); + Bfree(bs); + Bfree(bd0); + Bfree(delta); +ret: + if (se) + *se = (char *)s; + result = sign ? -value(rv) : value(rv); + + Bigint *&p5s = s_bigint_data->p5s; + while (p5s) { + tmp = p5s; + p5s = p5s->next; + free(tmp); + } + + return result; +} + double zend_hex_strtod(const char *str, char **endptr) { const char *s = str; diff --git a/hphp/runtime/base/zend/zend_strtod.h b/hphp/runtime/base/zend/zend_strtod.h index fc1f0fc4a..e0a86af9d 100644 --- a/hphp/runtime/base/zend/zend_strtod.h +++ b/hphp/runtime/base/zend/zend_strtod.h @@ -26,6 +26,7 @@ namespace HPHP { void zend_freedtoa(char *s); char * zend_dtoa(double _d, int mode, int ndigits, int *decpt, int *sign, char **rve); +double zend_strtod(const char *s00, char **se); double zend_hex_strtod(const char *str, char **endptr); double zend_oct_strtod(const char *str, char **endptr); int zend_startup_strtod(void);