From 202ac6aaf3f7e1e4adf0755647e714bd5f6df25c Mon Sep 17 00:00:00 2001
From: Sam Hocevar <sam@hocevar.net>
Date: Wed, 21 Sep 2011 17:10:21 +0000
Subject: [PATCH] core: switch real mantissa to uint16_t instead of uint32_t to
 ease the multiplication.

---
 src/real.cpp       | 74 +++++++++++++++++++++++++++++++++++++++++-----
 src/real.h         | 16 ++++++----
 test/unit/real.cpp | 24 +++++++++++++--
 3 files changed, 99 insertions(+), 15 deletions(-)

diff --git a/src/real.cpp b/src/real.cpp
index dfe0bd52..ee5fee50 100644
--- a/src/real.cpp
+++ b/src/real.cpp
@@ -13,6 +13,7 @@
 #endif
 
 #include <cstring>
+#include <cstdio>
 
 #include "core.h"
 
@@ -21,28 +22,87 @@ using namespace std;
 namespace lol
 {
 
-template<> real4k::Real(float f)
+real::real(float f)
 {
     union { float f; uint32_t x; } u = { f };
 
     uint32_t sign = u.x & 0x80000000u;
-    int e = ((u.x >> 23) & 0xff) + (1 << 30) - (1 << 10);
+    int e = ((u.x >> 23) & 0xff) + (1 << 30) - (1 << 7);
 
     m_signexp = sign | e;
-    m_mantissa[0] = u.x << 17;
-    memset(m_mantissa + 1, 0, sizeof(m_mantissa) - sizeof(m_mantissa[0]));
+    m_mantissa[0] = u.x >> 7;
+    m_mantissa[1] = u.x << 9;
+    memset(m_mantissa + 2, 0, sizeof(m_mantissa) - sizeof(m_mantissa[0]));
 }
 
-template<> real4k::operator float() const
+real::operator float() const
 {
     union { float f; uint32_t x; } u;
 
-    u.x = m_mantissa[0] >> 17;
-    u.x |= ((m_signexp & 0x7fffffffu) - (1 << 30) + (1 << 10)) << 23;
+    u.x = m_mantissa[0] << 7;
+    u.x |= m_mantissa[1] >> 9;
+    u.x |= ((m_signexp & 0x7fffffffu) - (1 << 30) + (1 << 7)) << 23;
     u.x |= m_signexp & 0x80000000u;
 
     return u.f;
 }
 
+real real::operator *(real const &x) const
+{
+    real ret;
+
+    ret.m_signexp = (m_signexp ^ x.m_signexp) & 0x80000000u;
+    int e = (m_signexp & 0x7fffffffu) - (1 << 30) + 1
+          + (x.m_signexp & 0x7fffffffu) - (1 << 30) + 1;
+
+    /* Accumulate low order product; no need to store it, we just
+     * want the carry value */
+    uint32_t carry = 0;
+    for (int i = 0; i < BIGITS; i++)
+    {
+        for (int j = 0; j < i + 1; j++)
+            carry += m_mantissa[BIGITS - 1 - j]
+                   * x.m_mantissa[BIGITS - 1 + j - i];
+        carry >>= 16;
+    }
+
+    for (int i = 0; i < BIGITS; i++)
+    {
+        for (int j = i + 1; j < BIGITS; j++)
+            carry += m_mantissa[BIGITS - 1 - j]
+                   * x.m_mantissa[j - 1 - i];
+
+        carry += m_mantissa[BIGITS - 1 - i];
+        carry += x.m_mantissa[BIGITS - 1 - i];
+        ret.m_mantissa[BIGITS - 1 - i] = carry & 0xffffu;
+        carry >>= 16;
+    }
+
+    /* Renormalise in case we overflowed the mantissa */
+    if (carry)
+    {
+        carry--;
+        for (int i = 0; i < BIGITS; i++)
+        {
+            uint16_t tmp = ret.m_mantissa[i];
+            ret.m_mantissa[i] = (carry << 15) | (tmp >> 1);
+            carry = tmp & 0x0001u;
+        }
+        e++;
+    }
+
+    ret.m_signexp |= e + (1 << 30) - 1;
+
+    return ret;
+}
+
+void real::print() const
+{
+    printf("%x  %08x  ", m_signexp >> 31, (m_signexp << 1) >> 1);
+    for (int i = 0; i < BIGITS; i++)
+        printf("%04x ", m_mantissa[i]);
+    printf("\n");
+}
+
 } /* namespace lol */
 
diff --git a/src/real.h b/src/real.h
index d31238c3..15b7a45e 100644
--- a/src/real.h
+++ b/src/real.h
@@ -21,21 +21,25 @@
 namespace lol
 {
 
-template<int NBITS> class Real
+class real
 {
+    static int const BIGITS = 10;
+
 public:
-    inline Real<NBITS>() {}
-    Real<NBITS>(float f);
+    inline real() { }
+    real(float f);
 
     operator float() const;
+    real operator *(real const &x) const;
+
+    void print() const;
 
 private:
+    uint32_t m_size;
     uint32_t m_signexp;
-    uint32_t m_mantissa[(NBITS + 31) / 32];
+    uint16_t m_mantissa[BIGITS];
 };
 
-typedef Real<4096> real4k;
-
 } /* namespace lol */
 
 #endif // __LOL_REAL_H__
diff --git a/test/unit/real.cpp b/test/unit/real.cpp
index 18527e1d..6af9122b 100644
--- a/test/unit/real.cpp
+++ b/test/unit/real.cpp
@@ -25,11 +25,31 @@ LOLUNIT_FIXTURE(RealTest)
 public:
     LOLUNIT_TEST(test_real_from_float)
     {
-        float x = real4k(0.0f);
-        float y = real4k(1.0f);
+        float x = real(0.0f);
+        float y = real(1.0f);
+        float z = real(1.5f);
 
         LOLUNIT_ASSERT_EQUAL(x, 0.0f);
         LOLUNIT_ASSERT_EQUAL(y, 1.0f);
+        LOLUNIT_ASSERT_EQUAL(z, 1.5f);
+    }
+
+    LOLUNIT_TEST(test_real_mul)
+    {
+        real x(1.25f);
+        real y(1.5f);
+        real z(1.99999f);
+        real w(-1.5f);
+
+        float m1 = x * x;
+        float m2 = y * y;
+        float m3 = z * z;
+        float m4 = w * w;
+
+        LOLUNIT_ASSERT_EQUAL(m1, 1.25f * 1.25f);
+        LOLUNIT_ASSERT_EQUAL(m2, 1.5f * 1.5f);
+        LOLUNIT_ASSERT_EQUAL(m3, 1.99999f * 1.99999f);
+        LOLUNIT_ASSERT_EQUAL(m4, -1.5f * -1.5f);
     }
 };