From faa877ff7818925995e682e40b3ed491c717da83 Mon Sep 17 00:00:00 2001 From: Milo Yip Date: Fri, 19 Sep 2014 08:59:36 +0800 Subject: [PATCH] Partial StrtodDiyFp implementation [ci skip] --- include/rapidjson/internal/diyfp.h | 25 +++++- include/rapidjson/internal/strtod.h | 113 ++++++++++++++++++++-------- test/unittest/readertest.cpp | 2 +- 3 files changed, 106 insertions(+), 34 deletions(-) diff --git a/include/rapidjson/internal/diyfp.h b/include/rapidjson/internal/diyfp.h index afd62bab..19aedb19 100644 --- a/include/rapidjson/internal/diyfp.h +++ b/include/rapidjson/internal/diyfp.h @@ -155,7 +155,7 @@ struct DiyFp { int e; }; -inline DiyFp GetCachedPower(int e, int* K) { +inline uint64_t GetCachedPowerF(size_t index) { // 10^-348, 10^-340, ..., 10^340 static const uint64_t kCachedPowers_F[] = { RAPIDJSON_UINT64_C2(0xfa8fd5a0, 0x081c0288), RAPIDJSON_UINT64_C2(0xbaaee17f, 0xa23ebf76), @@ -203,6 +203,10 @@ inline DiyFp GetCachedPower(int e, int* K) { RAPIDJSON_UINT64_C2(0x9e19db92, 0xb4e31ba9), RAPIDJSON_UINT64_C2(0xeb96bf6e, 0xbadf77d9), RAPIDJSON_UINT64_C2(0xaf87023b, 0x9bf0ee6b) }; + return kCachedPowers_F[index]; +} + +inline DiyFp GetCachedPower(int e, int* K) { static const int16_t kCachedPowers_E[] = { -1220, -1193, -1166, -1140, -1113, -1087, -1060, -1034, -1007, -980, -954, -927, -901, -874, -847, -821, -794, -768, -741, -715, @@ -224,9 +228,26 @@ inline DiyFp GetCachedPower(int e, int* K) { unsigned index = static_cast((k >> 3) + 1); *K = -(-348 + static_cast(index << 3)); // decimal exponent no need lookup table - return DiyFp(kCachedPowers_F[index], kCachedPowers_E[index]); + return DiyFp(GetCachedPowerF(index), kCachedPowers_E[index]); } +inline DiyFp GetCachedPower10(int exp, int *outExp) { + // static const int16_t kCachedPowers_E10[] = { + // -348, -340, -332, -324, -316, -308, -300, -292, -284, -276, + // -268, -260, -252, -244, -236, -228, -220, -212, -204, -196, + // -188, -180, -172, -164, -156, -148, -140, -132, -124, -116, + // -108, -100, -92, -84, -76, -68, -60, -52, -44, -36, + // -28, -20, -12, -4, 4, 12, 20, 28, 36, 44, + // 52, 60, 68, 76, 84, 92, 100, 108, 116, 124, + // 132, 140, 148, 156, 164, 172, 180, 188, 196, 204, + // 212, 220, 228, 236, 244, 252, 260, 268, 276, 284, + // 292, 300, 308, 316, 324, 332, 340 + // }; + unsigned index = (exp + 348) / 8; + *outExp = -348 + index * 8; + return GetCachedPowerF(index); + } + #ifdef __GNUC__ RAPIDJSON_DIAG_POP #endif diff --git a/include/rapidjson/internal/strtod.h b/include/rapidjson/internal/strtod.h index e1c09e17..b18f19e2 100644 --- a/include/rapidjson/internal/strtod.h +++ b/include/rapidjson/internal/strtod.h @@ -24,6 +24,7 @@ #include "../rapidjson.h" #include "ieee754.h" #include "biginteger.h" +#include "diyfp.h" #include "pow10.h" namespace rapidjson { @@ -141,7 +142,84 @@ inline bool StrtodFast(double d, int p, double* result) { return false; } -inline double StrtodBigInteger(double d, int p, const char* decimals, size_t length, size_t decimalPosition, int exp) { +// Compute an approximation and see if it is within 1/2 ULP +inline bool StrtodDiyFp(const char* decimals, size_t length, size_t decimalPosition, int exp, double* result) { + uint64_t significand = 0; + size_t i = 0; // 2^64 - 1 = 18446744073709551615, 1844674407370955161 = 0x199999990x99999999 + for (; i < length && (significand <= RAPIDJSON_UINT64_C2(0x19999999, 0x99999999) || decimals[i] <= '4'); i++) + significand = significand * 10 + (decimals[i] - '0'); + + if (i < length && decimals[i] >= '5') // Rounding + significand++; + + DiyFp v(significand, 0); + size_t remaining = length - i; + const int dExp = (int)decimalPosition - i + exp; + exp += (int)remaining; + + const unsigned kUlpShift = 3; + const unsigned kUlp = 1 << kUlpShift; + int error = (remaining == 0) ? 0 : kUlp / 2; + + v = v.Normalize(); + error <<= - v.e; + + int actualExp; + v = v * GetCachedPower10(exp, &actualExp); + if (actualExp != exp) { + static const DiyFp kPow10[] = { + DiyFp(RAPIDJSON_UINT64_C2(0xa0000000, 00000000), -60), // 10^1 + DiyFp(RAPIDJSON_UINT64_C2(0xc8000000, 00000000), -57), // 10^2 + DiyFp(RAPIDJSON_UINT64_C2(0xfa000000, 00000000), -54), // 10^3 + DiyFp(RAPIDJSON_UINT64_C2(0x9c400000, 00000000), -50), // 10^4 + DiyFp(RAPIDJSON_UINT64_C2(0xc3500000, 00000000), -47), // 10^5 + DiyFp(RAPIDJSON_UINT64_C2(0xf4240000, 00000000), -44), // 10^6 + DiyFp(RAPIDJSON_UINT64_C2(0x98968000, 00000000), -40) // 10^7 + }; + v = v * kPow10[actualExp - exp - 1]; + } + + error += kUlp + (error == 0 ? 0 : 1); + + int oldExp = v.e; + v = v.Normalize(); + error <<= oldExp - v.e; + + return true; +} + +inline double StrtodBigInteger(double approx, const char* decimals, size_t length, size_t decimalPosition, int exp) { + const BigInteger dInt(decimals, length); + const int dExp = (int)decimalPosition - (int)length + exp; + Double a(approx); + for (int i = 0; i < 10; i++) { + bool adjustToNegative; + int cmp = CheckWithinHalfULP(a.Value(), dInt, dExp, &adjustToNegative); + if (cmp < 0) + return a.Value(); // within half ULP + else if (cmp == 0) { + // Round towards even + if (a.Significand() & 1) + return adjustToNegative ? a.PreviousPositiveDouble() : a.NextPositiveDouble(); + else + return a.Value(); + } + else // adjustment + a = adjustToNegative ? a.PreviousPositiveDouble() : a.NextPositiveDouble(); + } + + // This should not happen, but in case there is really a bug, break the infinite-loop + return a.Value(); +} + +inline double StrtodFullPrecision(double d, int p, const char* decimals, size_t length, size_t decimalPosition, int exp) { + RAPIDJSON_ASSERT(d >= 0.0); + RAPIDJSON_ASSERT(length >= 1); + + double result; + if (StrtodFast(d, p, &result)) + return result; + // Trim leading zeros while (*decimals == '0' && length > 1) { length--; @@ -167,38 +245,11 @@ inline double StrtodBigInteger(double d, int p, const char* decimals, size_t len if (int(length) + exp < -324) return 0.0; - const BigInteger dInt(decimals, length); - const int dExp = (int)decimalPosition - (int)length + exp; - Double approx = StrtodNormalPrecision(d, p); - for (int i = 0; i < 10; i++) { - bool adjustToNegative; - int cmp = CheckWithinHalfULP(approx.Value(), dInt, dExp, &adjustToNegative); - if (cmp < 0) - return approx.Value(); // within half ULP - else if (cmp == 0) { - // Round towards even - if (approx.Significand() & 1) - return adjustToNegative ? approx.PreviousPositiveDouble() : approx.NextPositiveDouble(); - else - return approx.Value(); - } - else // adjustment - approx = adjustToNegative ? approx.PreviousPositiveDouble() : approx.NextPositiveDouble(); - } - - // This should not happen, but in case there is really a bug, break the infinite-loop - return approx.Value(); -} - -inline double StrtodFullPrecision(double d, int p, const char* decimals, size_t length, size_t decimalPosition, int exp) { - RAPIDJSON_ASSERT(d >= 0.0); - RAPIDJSON_ASSERT(length >= 1); - - double result; - if (StrtodFast(d, p, &result)) + if (StrtodDiyFp(decimals, length, decimalPosition, exp, &result)) return result; - return StrtodBigInteger(d, p, decimals, length, decimalPosition, exp); + // Use approximation from StrtodDiyFp and make adjustment with BigInteger comparison + return StrtodBigInteger(result, decimals, length, decimalPosition, exp); } } // namespace internal diff --git a/test/unittest/readertest.cpp b/test/unittest/readertest.cpp index 7a449567..882327a3 100644 --- a/test/unittest/readertest.cpp +++ b/test/unittest/readertest.cpp @@ -258,7 +258,7 @@ static void TestParseDouble() { } #if 1 - static const unsigned count = 1000000; + static const unsigned count = 10000000; // Random test for double { Random r;