From 93479c28760e7cb980a49e72b994ac1babc306ac Mon Sep 17 00:00:00 2001
From: Sam Hocevar <sam@hocevar.net>
Date: Thu, 22 Sep 2011 07:16:05 +0000
Subject: [PATCH] core: implement real subtraction.

---
 src/real.cpp       | 78 +++++++++++++++++++++++++++++++++++++++++++++-
 test/unit/real.cpp | 13 ++------
 2 files changed, 80 insertions(+), 11 deletions(-)

diff --git a/src/real.cpp b/src/real.cpp
index 548b409a..544bfa92 100644
--- a/src/real.cpp
+++ b/src/real.cpp
@@ -103,7 +103,7 @@ real real::operator +(real const &x) const
     uint32_t carry = 0;
     for (int i = BIGITS; i--; )
     {
-        carry = m_mantissa[i];
+        carry += m_mantissa[i];
         if (i - bigoff >= 0)
             carry += x.m_mantissa[i - bigoff] >> off;
         else if (i - bigoff == -1)
@@ -153,6 +153,82 @@ real real::operator -(real const &x) const
 
     real ret;
 
+    int e1 = m_signexp - (1 << 30) + 1;
+    int e2 = x.m_signexp - (1 << 30) + 1;
+
+    int bigoff = (e1 - e2) / (sizeof(uint16_t) * 8);
+    int off = e1 - e2 - bigoff * (sizeof(uint16_t) * 8);
+
+    ret.m_signexp = m_signexp;
+
+    int32_t carry = 0;
+    for (int i = 0; i < bigoff; i++)
+    {
+        carry -= x.m_mantissa[BIGITS - i];
+        carry = (carry & 0xffff0000u) | (carry >> 16);
+    }
+    carry -= x.m_mantissa[BIGITS - 1 - bigoff] & ((1 << off) - 1);
+    carry /= (1 << off);
+
+    for (int i = BIGITS; i--; )
+    {
+        carry += m_mantissa[i];
+        if (i - bigoff >= 0)
+            carry -= x.m_mantissa[i - bigoff] >> off;
+        else if (i - bigoff == -1)
+            carry -= 0x0001u >> off;
+
+        if (i - bigoff > 0)
+            carry -= (x.m_mantissa[i - bigoff - 1] << (16 - off)) & 0xffffu;
+        else if (i - bigoff == 0)
+            carry -= 0x0001u << (16 - off);
+
+        ret.m_mantissa[i] = carry;
+        carry = (carry & 0xffff0000u) | (carry >> 16);
+    }
+
+    carry += 1;
+
+    /* Renormalise if we underflowed the mantissa */
+    if (carry == 0)
+    {
+        /* How much do we need to shift the mantissa? FIXME: this could
+         * be computed above */
+        off = 0;
+        for (int i = 0; i < BIGITS; i++)
+        {
+            if (!ret.m_mantissa[i])
+            {
+                off += sizeof(uint16_t) * 8;
+                continue;
+            }
+
+            for (uint16_t tmp = ret.m_mantissa[i]; tmp < 0x8000u; tmp <<= 1)
+                off++;
+            break;
+        }
+        if (off == BIGITS * sizeof(uint16_t) * 8)
+            ret.m_signexp &= 0x80000000u;
+        else
+        {
+            off++; /* Shift one more to get rid of the leading one */
+            ret.m_signexp -= off;
+
+            bigoff = off / (sizeof(uint16_t) * 8);
+            off -= bigoff * sizeof(uint16_t) * 8;
+
+            for (int i = 0; i < BIGITS; i++)
+            {
+                uint16_t tmp = 0;
+                if (i + bigoff < BIGITS)
+                    tmp |= ret.m_mantissa[i + bigoff] << off;
+                if (i + bigoff + 1 < BIGITS)
+                    tmp |= ret.m_mantissa[i + bigoff + 1] >> (16 - off);
+                ret.m_mantissa[i] = tmp;
+            }
+        }
+    }
+
     return ret;
 }
 
diff --git a/test/unit/real.cpp b/test/unit/real.cpp
index d84e49c7..816815bf 100644
--- a/test/unit/real.cpp
+++ b/test/unit/real.cpp
@@ -88,16 +88,9 @@ LOLUNIT_FIXTURE(RealTest)
 
     LOLUNIT_TEST(test_real_sub)
     {
-#if 0
-printf("\n");
-real k(1.25f);
-k.print();
-real l(1.0f);
-l.print();
-real m = k - l;
-m.print();
-#endif
-        LOLUNIT_ASSERT(true);
+        float a1 = real(1.0f) + real(1e20f) - real(1e20f);
+
+        LOLUNIT_ASSERT_EQUAL(a1, 1.0f);
     }
 
     LOLUNIT_TEST(test_real_mul)