Changeset 975Tweet

Timestamp:

Sep 22, 2011, 6:16:42 PM (11 years ago)

Author:

sam

Message:

core: implement division of reals and change their default precision
to 32 bigits.

Location:

trunk

Files:

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

trunk/src/real.cpp

@@ -282 +282 @@
 }
 
+real real::operator /(real const &x) const
+{
+    return *this * fres(x);
+}
+
 bool real::operator <(real const &x) const
 {
@@ -332 +337 @@
 }
 
+real fres(real const &x)
+{
+    if (!(x.m_signexp << 1))
+    {
+        real ret = x;
+        ret.m_signexp = x.m_signexp | 0x7fffffffu;
+        ret.m_mantissa[0] = 0;
+        return ret;
+    }
+
+    /* Use the system's float inversion to approximate 1/x */
+    union { float f; uint32_t x; } u = { 1.0f }, v = { 1.0f };
+    v.x |= (uint32_t)x.m_mantissa[0] << 7;
+    v.x |= (uint32_t)x.m_mantissa[1] >> 9;
+    v.f = 1.0 / v.f;
+
+    real ret;
+    ret.m_mantissa[0] = (v.x >> 7) & 0xffffu;
+    ret.m_mantissa[1] = (v.x << 9) & 0xffffu;
+    /* Better convergence with the mantissa zeroed. */
+    memset(ret.m_mantissa + 2, 0, (real::BIGITS - 2) * sizeof(uint16_t));
+
+    uint32_t sign = x.m_signexp & 0x80000000u;
+    ret.m_signexp = sign;
+
+    int exponent = (x.m_signexp & 0x7fffffffu) + 1;
+    exponent = -exponent + (v.x >> 23) - (u.x >> 23);
+    ret.m_signexp |= (exponent - 1) & 0x7fffffffu;
+
+    /* Five steps of Newton-Raphson seems enough for 32-bigit reals. */
+    real two(2.0f);
+    ret = ret * (two - ret * x);
+    ret = ret * (two - ret * x);
+    ret = ret * (two - ret * x);
+    ret = ret * (two - ret * x);
+    ret = ret * (two - ret * x);
+
+    return ret;
+}
+
 void real::print() const
 {

trunk/src/real.h

@@ -24 +24 @@
 class real
 {
-    static int const BIGITS = 10;
+    static int const BIGITS = 32;
 
 public:
@@ -31 +31 @@
 
     operator float() const;
+
     real operator -() const;
     real operator +(real const &x) const;
     real operator -(real const &x) const;
     real operator *(real const &x) const;
+    real operator /(real const &x) const;
 
     bool operator <(real const &x) const;
@@ -40 +42 @@
     bool operator <=(real const &x) const;
     bool operator >=(real const &x) const;
+
+    friend real fres(real const &x);
 
     void print() const;

trunk/test/unit/real.cpp

@@ -111 +111 @@
         LOLUNIT_ASSERT_EQUAL(m4, -1.5f * -1.5f);
     }
+
+    LOLUNIT_TEST(test_real_div)
+    {
+        real a1(1.0f);
+        real a2(2.0f);
+
+        float m1 = a1 / a1;
+        float m2 = a2 / a1;
+        float m3 = a1 / a2;
+        float m4 = a2 / a2;
+
+        LOLUNIT_ASSERT_EQUAL(m1, 1.0f);
+        LOLUNIT_ASSERT_EQUAL(m2, 2.0f);
+        LOLUNIT_ASSERT_EQUAL(m3, 0.5f);
+        LOLUNIT_ASSERT_EQUAL(m4, 1.0f);
+    }
 };