From 5d9167bda05e12c6defd5d1958f278277049bea0 Mon Sep 17 00:00:00 2001
From: Sam Hocevar <sam@hocevar.net>
Date: Tue, 11 Oct 2011 22:01:44 +0000
Subject: [PATCH] core: fix an accuracy error in real::re() and real::sqrt()
 introduced in the 16-to-32-bit refactoring.

---
 src/real.cpp       | 16 +++++++---------
 src/real.h         |  3 ++-
 test/unit/real.cpp | 25 ++++++++++++++++---------
 3 files changed, 25 insertions(+), 19 deletions(-)

diff --git a/src/real.cpp b/src/real.cpp
index fbb632e2..b6427317 100644
--- a/src/real.cpp
+++ b/src/real.cpp
@@ -22,8 +22,6 @@ using namespace std;
 namespace lol
 {
 
-static int const BIGIT_BITS = 32;
-
 real::real(float f) { *this = (double)f; }
 real::real(int i) { *this = (double)i; }
 real::real(unsigned int i) { *this = (double)i; }
@@ -500,9 +498,9 @@ real re(real const &x)
     exponent = -exponent + (v.x >> 23) - (u.x >> 23);
     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);
 
     return ret;
@@ -544,9 +542,9 @@ real sqrt(real const &x)
     exponent = exponent + (v.x >> 23) - (u.x >> 23);
     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);
         ret.m_signexp--;
@@ -668,10 +666,10 @@ real floor(real const &x)
     {
         if (exponent <= 0)
             ret.m_mantissa[i] = 0;
-        else if (exponent < BIGIT_BITS)
-            ret.m_mantissa[i] &= ~((1 << (BIGIT_BITS - exponent)) - 1);
+        else if (exponent < real::BIGIT_BITS)
+            ret.m_mantissa[i] &= ~((1 << (real::BIGIT_BITS - exponent)) - 1);
 
-        exponent -= BIGIT_BITS;
+        exponent -= real::BIGIT_BITS;
     }
 
     return ret;
diff --git a/src/real.h b/src/real.h
index ab331f28..e4fab283 100644
--- a/src/real.h
+++ b/src/real.h
@@ -104,12 +104,13 @@ public:
     static real const R_SQRT3;
     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;
     uint32_t m_mantissa[BIGITS];
diff --git a/test/unit/real.cpp b/test/unit/real.cpp
index 93e17120..323ce4f6 100644
--- a/test/unit/real.cpp
+++ b/test/unit/real.cpp
@@ -180,16 +180,13 @@ LOLUNIT_FIXTURE(RealTest)
         LOLUNIT_ASSERT_EQUAL(m4, -1.5f * -1.5f);
     }
 
-    LOLUNIT_TEST(Division)
+    LOLUNIT_TEST(ExactDivision)
     {
-        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;
+        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);
         LOLUNIT_ASSERT_EQUAL(m2, 2.0f);
@@ -198,6 +195,16 @@ LOLUNIT_FIXTURE(RealTest)
         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);
+    }
+
     LOLUNIT_TEST(Shift)
     {
         real a1(1.5);