Changeset 1017Tweet

Timestamp:

Oct 12, 2011, 12:01:44 AM (11 years ago)

Author:

sam

Message:

core: fix an accuracy error in real::re() and real::sqrt() introduced in
the 16-to-32-bit refactoring.

Location:

trunk

Files:

• src/real.cpp (modified) (4 diffs)
• src/real.h (modified) (1 diff)
• test/unit/real.cpp (modified) (2 diffs)

trunk/src/real.cpp

@@ -23 +23 @@
 {
 
-static int const BIGIT_BITS = 32;
-
 real::real(float f) { *this = (double)f; }
 real::real(int i) { *this = (double)i; }
@@ -501 +499 @@
     ret.m_signexp |= (exponent - 1) & 0x7fffffffu;
 
-    /* FIXME: log2(BIGITS) steps of Newton-Raphson seems to be enough for
+    /* FIXME: 1+log2(BIGITS) steps of Newton-Raphson seems to be enough for
      * convergence, but this hasn't been checked seriously. */
-    for (int i = 2; i < real::BIGITS; i *= 2)
+    for (int i = 1; i <= real::BIGITS; i *= 2)
         ret = ret * (real::R_2 - ret * x);
 
@@ -545 +543 @@
     ret.m_signexp |= exponent & 0x7fffffffu;
 
-    /* FIXME: log2(BIGITS) steps of Newton-Raphson seems to be enough for
+    /* FIXME: 1+log2(BIGITS) steps of Newton-Raphson seems to be enough for
      * convergence, but this hasn't been checked seriously. */
-    for (int i = 2; i < real::BIGITS; i *= 2)
+    for (int i = 1; i <= real::BIGITS; i *= 2)
     {
         ret = ret * (real::R_3 - ret * ret * x);
@@ -669 +667 @@
         if (exponent <= 0)
             ret.m_mantissa[i] = 0;
-        else if (exponent < BIGIT_BITS)
-            ret.m_mantissa[i] &= ~((1 << (BIGIT_BITS - exponent)) - 1);
-
-        exponent -= BIGIT_BITS;
+        else if (exponent < real::BIGIT_BITS)
+            ret.m_mantissa[i] &= ~((1 << (real::BIGIT_BITS - exponent)) - 1);
+
+        exponent -= real::BIGIT_BITS;
     }
 

trunk/src/real.h

@@ -105 +105 @@
     static real const R_SQRT1_2;
 
-private:
     /* XXX: changing this requires tuning real::fres (the number of
      * Newton-Raphson iterations) and real::print (the number of printed
      * digits) */
     static int const BIGITS = 16;
+    static int const BIGIT_BITS = 32;
 
+private:
     uint32_t m_size;
     uint32_t m_signexp;

trunk/test/unit/real.cpp

@@ -181 +181 @@
     }
 
-    LOLUNIT_TEST(Division)
-    {
-        real a1(1.0f);
-        real a2(2.0f);
-
-        float m1 = a1 / a1;
-        float m2 = a2 / a1;
-        float m3 = a1 / a2;
-        float m4 = a2 / a2;
-        float m5 = a1 / -a2;
+    LOLUNIT_TEST(ExactDivision)
+    {
+        float m1 = real::R_1 / real::R_1;
+        float m2 = real::R_2 / real::R_1;
+        float m3 = real::R_1 / real::R_2;
+        float m4 = real::R_2 / real::R_2;
+        float m5 = real::R_1 / -real::R_2;
 
         LOLUNIT_ASSERT_EQUAL(m1, 1.0f);
@@ -197 +194 @@
         LOLUNIT_ASSERT_EQUAL(m4, 1.0f);
         LOLUNIT_ASSERT_EQUAL(m5, -0.5f);
+    }
+
+    LOLUNIT_TEST(InexactDivision)
+    {
+        /* 1 / 3 * 3 should be close to 1... check that it does not differ
+         * by more than 2^-k where k is the number of bits in the mantissa. */
+        real a = real::R_1 / real::R_3 * real::R_3;
+        real b = (real::R_1 - a) << (real::BIGITS * real::BIGIT_BITS);
+
+        LOLUNIT_ASSERT_LEQUAL((double)fabs(b), 1.0);
     }