Partial StrtodDiyFp implementation [ci skip]

2 files changed, 104 insertions, 32 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<unsigned>((k >> 3) + 1);
     *K = -(-348 + static_cast<int>(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,53 +142,74 @@ 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) {
-    // Trim leading zeros
-    while (*decimals == '0' && length > 1) {
-        length--;
-        decimals++;
-        decimalPosition--;
-    }
+// 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++;
 
-    // Trim trailing zeros
-    while (decimals[length - 1] == '0' && length > 1) {
-        length--;
-        decimalPosition--;
-        exp++;
-    }
+    DiyFp v(significand, 0);
+    size_t remaining = length - i;
+    const int dExp = (int)decimalPosition - i + exp;
+    exp += (int)remaining;
 
-    // Trim right-most digits
-    const int kMaxDecimalDigit = 780;
-    if (length > kMaxDecimalDigit) {
-        exp += (int(length) - kMaxDecimalDigit);
-        length = kMaxDecimalDigit;
+    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];
     }
 
-    // If too small, underflow to zero
-    if (int(length) + exp < -324)
-        return 0.0;
+    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 approx = StrtodNormalPrecision(d, p);
+    Double a(approx);
     for (int i = 0; i < 10; i++) {
         bool adjustToNegative;
-        int cmp = CheckWithinHalfULP(approx.Value(), dInt, dExp, &adjustToNegative);
+        int cmp = CheckWithinHalfULP(a.Value(), dInt, dExp, &adjustToNegative);
         if (cmp < 0)
-            return approx.Value();  // within half ULP
+            return a.Value();  // within half ULP
         else if (cmp == 0) {
             // Round towards even
-            if (approx.Significand() & 1)
-                return adjustToNegative ? approx.PreviousPositiveDouble() : approx.NextPositiveDouble();
+            if (a.Significand() & 1)
+                return adjustToNegative ? a.PreviousPositiveDouble() : a.NextPositiveDouble();
             else
-                return approx.Value();
+                return a.Value();
         }
         else // adjustment
-            approx = adjustToNegative ? approx.PreviousPositiveDouble() : approx.NextPositiveDouble();
+            a = adjustToNegative ? a.PreviousPositiveDouble() : a.NextPositiveDouble();
     }
 
     // This should not happen, but in case there is really a bug, break the infinite-loop
-    return approx.Value();
+    return a.Value();
 }
 
 inline double StrtodFullPrecision(double d, int p, const char* decimals, size_t length, size_t decimalPosition, int exp) {
@@ -198,7 +220,36 @@ inline double StrtodFullPrecision(double d, int p, const char* decimals, size_t
     if (StrtodFast(d, p, &result))
         return result;
 
-    return StrtodBigInteger(d, p, decimals, length, decimalPosition, exp);
+    // Trim leading zeros
+    while (*decimals == '0' && length > 1) {
+        length--;
+        decimals++;
+        decimalPosition--;
+    }
+
+    // Trim trailing zeros
+    while (decimals[length - 1] == '0' && length > 1) {
+        length--;
+        decimalPosition--;
+        exp++;
+    }
+
+    // Trim right-most digits
+    const int kMaxDecimalDigit = 780;
+    if (length > kMaxDecimalDigit) {
+        exp += (int(length) - kMaxDecimalDigit);
+        length = kMaxDecimalDigit;
+    }
+
+    // If too small, underflow to zero
+    if (int(length) + exp < -324)
+        return 0.0;
+
+    if (StrtodDiyFp(decimals, length, decimalPosition, exp, &result))
+        return result;
+
+    // Use approximation from StrtodDiyFp and make adjustment with BigInteger comparison
+    return StrtodBigInteger(result, decimals, length, decimalPosition, exp);
 }
 
 } // namespace internal