Commit cbd74752 authored by Milo Yip's avatar Milo Yip

Fix normal-subnormal boundary and add more boundary cases in unit tests.

parent fa52a269
......@@ -122,6 +122,22 @@ public:
return *this;
}
BigInteger& operator+=(const BigInteger& rhs) {
size_t count = count_ > rhs.count_ ? count_ : rhs.count_;
bool carry = false;
for (size_t i = 0; i < count; i++) {
bool outCarry;
digits_[i] = FullAdd(i < count_ ? digits_[i] : 0, i < rhs.count_ ? rhs.digits_[i] : 0, carry, &outCarry);
carry = outCarry;
}
count_ = count;
if (carry)
PushBack(1);
return *this;
}
BigInteger& operator*=(uint64_t u) {
if (u == 0) return *this = 0;
if (u == 1) return *this;
......@@ -332,6 +348,12 @@ private:
#endif
}
static Type FullAdd(Type a, Type b, bool inCarry, bool* outCarry) {
Type c = a + b + (inCarry ? 1 : 0);
*outCarry = c < a;
return c;
}
static const size_t kBitCount = 3328; // 64bit * 54 > 10^1000
static const size_t kCapacity = kBitCount / sizeof(Type);
static const size_t kTypeBit = sizeof(Type) * 8;
......@@ -429,16 +451,14 @@ inline int CheckWithinHalfULP(double b, const BigInteger& d, int dExp, bool* adj
*adjustToNegative = dS.Difference(bS, &delta);
int cmp = delta.Compare(hS);
// If delta is within 1/2 ULP, check for special case when significand is power of two.
// In this case, need to compare with 1/2h in the lower bound.
if (cmp < 0 && *adjustToNegative && db.IsNormal() && (bInt & (bInt - 1)) == 0) {
if (cmp < 0 && *adjustToNegative && // within and dS < bS
db.IsNormal() && (bInt & (bInt - 1)) == 0 && // Power of 2
db.Uint64Value() != RAPIDJSON_UINT64_C2(0x00100000, 0x00000000)) // minimum normal number must not do this
{
delta <<= 1;
if (delta.Compare(hS) <= 0)
return -1;
else
return 1;
return delta.Compare(hS);
}
return cmp;
}
......@@ -465,7 +485,6 @@ inline double FullPrecision(double d, int dExp, const char* decimals, size_t len
const BigInteger dInt(decimals, length);
Double approx = NormalPrecision(d, p);
bool lastAdjustToNegative = false;
for (int i = 0; i < 10; i++) {
bool adjustToNegative;
int cmp = CheckWithinHalfULP(approx.Value(), dInt, dExp, &adjustToNegative);
......@@ -479,17 +498,12 @@ inline double FullPrecision(double d, int dExp, const char* decimals, size_t len
return approx.Value();
}
else {
// If oscillate between negative/postive, terminate
if (i > 0 && adjustToNegative != lastAdjustToNegative)
break;
// adjustment
if (adjustToNegative)
approx = approx.PreviousDouble();
else
approx = approx.NextDouble();
}
lastAdjustToNegative = adjustToNegative;
}
// This should not happen, but in case there is really a bug, break the infinite-loop
......
......@@ -226,8 +226,27 @@ static void TestParseDouble() {
TEST_DOUBLE(fullPrecision, "123e34", 123e34); // Fast Path Cases In Disguise
TEST_DOUBLE(fullPrecision, "45913141877270640000.0", 45913141877270640000.0);
TEST_DOUBLE(fullPrecision, "2.2250738585072011e-308", 2.2250738585072011e-308); // http://www.exploringbinary.com/php-hangs-on-numeric-value-2-2250738585072011e-308/
//TEST_DOUBLE(fullPrecision, "2.2250738585072012e-308", 2.2250738585072012e-308); // http://www.exploringbinary.com/java-hangs-when-converting-2-2250738585072012e-308/
TEST_DOUBLE(fullPrecision, "4503599627370495.9999999999", 4503599627370495.9999999999); // slightly smaller than 2^52 for power of two case
// Since
// abs((2^-1022 - 2^-1074) - 2.2250738585072012e-308) = 3.109754131239141401123495768877590405345064751974375599... ¡Á 10^-324
// abs((2^-1022) - 2.2250738585072012e-308) = 1.830902327173324040642192159804623318305533274168872044... ¡Á 10 ^ -324
// So 2.2250738585072012e-308 should round to 2^-1022 = 2.2250738585072014e-308
TEST_DOUBLE(fullPrecision, "2.2250738585072012e-308", 2.2250738585072014e-308); // http://www.exploringbinary.com/java-hangs-when-converting-2-2250738585072012e-308/
// More closer to normal/subnormal boundary
// boundary = 2^-1022 - 2^-1075 = 2.225073858507201136057409796709131975934819546351645648... ¡Á 10^-308
TEST_DOUBLE(fullPrecision, "2.22507385850720113605740979670913197593481954635164564e-308", 2.2250738585072009e-308);
TEST_DOUBLE(fullPrecision, "2.22507385850720113605740979670913197593481954635164565e-308", 2.2250738585072014e-308);
// 1.0 is in (1.0 - 2^-54, 1.0 + 2^-53)
// 1.0 - 2^-54 = 0.999999999999999944488848768742172978818416595458984375
TEST_DOUBLE(fullPrecision, "0.999999999999999944488848768742172978818416595458984375", 1.0); // round to even
TEST_DOUBLE(fullPrecision, "0.999999999999999944488848768742172978818416595458984374", 0.99999999999999989); // previous double
TEST_DOUBLE(fullPrecision, "0.999999999999999944488848768742172978818416595458984376", 1.0); // next double
// 1.0 + 2^-53 = 1.00000000000000011102230246251565404236316680908203125
TEST_DOUBLE(fullPrecision, "1.00000000000000011102230246251565404236316680908203125", 1.0); // round to even
TEST_DOUBLE(fullPrecision, "1.00000000000000011102230246251565404236316680908203124", 1.0); // previous double
TEST_DOUBLE(fullPrecision, "1.00000000000000011102230246251565404236316680908203126", 1.00000000000000022); // next double
{
char n1e308[310]; // '1' followed by 308 '0'
......
......@@ -48,9 +48,6 @@ TEST(Strtod, BigInteger_Constructor) {
}
TEST(Strtod, BigInteger_AddUint64) {
const BigInteger kZero(0);
const BigInteger kOne(1);
BigInteger a = kZero;
a += 0u;
EXPECT_TRUE(kZero == a);
......@@ -69,6 +66,30 @@ TEST(Strtod, BigInteger_AddUint64) {
EXPECT_TRUE(BIGINTEGER_LITERAL("36893488147419103231") == b);
}
TEST(Strtod, BigInteger_Add) {
BigInteger a = kZero;
a += kZero;
EXPECT_TRUE(kZero == a);
a += kOne;
EXPECT_TRUE(kOne == a);
a += kOne;
EXPECT_TRUE(BigInteger(2) == a);
a = kUint64Max;
a += kOne;
EXPECT_TRUE(kTwo64 == a);
a = kOne;
a += kTwo64;
EXPECT_TRUE(BIGINTEGER_LITERAL("18446744073709551617") == a);
a = kTwo64;
a += kOne;
EXPECT_TRUE(BIGINTEGER_LITERAL("18446744073709551617") == a);
}
TEST(Strtod, BigInteger_MultiplyUint64) {
BigInteger a = kZero;
a *= static_cast <uint64_t>(0);
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment