Changeset 975


Ignore:
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:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/real.cpp

    r974 r975  
    282282}
    283283
     284real real::operator /(real const &x) const
     285{
     286    return *this * fres(x);
     287}
     288
    284289bool real::operator <(real const &x) const
    285290{
     
    332337}
    333338
     339real fres(real const &x)
     340{
     341    if (!(x.m_signexp << 1))
     342    {
     343        real ret = x;
     344        ret.m_signexp = x.m_signexp | 0x7fffffffu;
     345        ret.m_mantissa[0] = 0;
     346        return ret;
     347    }
     348
     349    /* Use the system's float inversion to approximate 1/x */
     350    union { float f; uint32_t x; } u = { 1.0f }, v = { 1.0f };
     351    v.x |= (uint32_t)x.m_mantissa[0] << 7;
     352    v.x |= (uint32_t)x.m_mantissa[1] >> 9;
     353    v.f = 1.0 / v.f;
     354
     355    real ret;
     356    ret.m_mantissa[0] = (v.x >> 7) & 0xffffu;
     357    ret.m_mantissa[1] = (v.x << 9) & 0xffffu;
     358    /* Better convergence with the mantissa zeroed. */
     359    memset(ret.m_mantissa + 2, 0, (real::BIGITS - 2) * sizeof(uint16_t));
     360
     361    uint32_t sign = x.m_signexp & 0x80000000u;
     362    ret.m_signexp = sign;
     363
     364    int exponent = (x.m_signexp & 0x7fffffffu) + 1;
     365    exponent = -exponent + (v.x >> 23) - (u.x >> 23);
     366    ret.m_signexp |= (exponent - 1) & 0x7fffffffu;
     367
     368    /* Five steps of Newton-Raphson seems enough for 32-bigit reals. */
     369    real two(2.0f);
     370    ret = ret * (two - ret * x);
     371    ret = ret * (two - ret * x);
     372    ret = ret * (two - ret * x);
     373    ret = ret * (two - ret * x);
     374    ret = ret * (two - ret * x);
     375
     376    return ret;
     377}
     378
    334379void real::print() const
    335380{
  • trunk/src/real.h

    r973 r975  
    2424class real
    2525{
    26     static int const BIGITS = 10;
     26    static int const BIGITS = 32;
    2727
    2828public:
     
    3131
    3232    operator float() const;
     33
    3334    real operator -() const;
    3435    real operator +(real const &x) const;
    3536    real operator -(real const &x) const;
    3637    real operator *(real const &x) const;
     38    real operator /(real const &x) const;
    3739
    3840    bool operator <(real const &x) const;
     
    4042    bool operator <=(real const &x) const;
    4143    bool operator >=(real const &x) const;
     44
     45    friend real fres(real const &x);
    4246
    4347    void print() const;
  • trunk/test/unit/real.cpp

    r974 r975  
    111111        LOLUNIT_ASSERT_EQUAL(m4, -1.5f * -1.5f);
    112112    }
     113
     114    LOLUNIT_TEST(test_real_div)
     115    {
     116        real a1(1.0f);
     117        real a2(2.0f);
     118
     119        float m1 = a1 / a1;
     120        float m2 = a2 / a1;
     121        float m3 = a1 / a2;
     122        float m4 = a2 / a2;
     123
     124        LOLUNIT_ASSERT_EQUAL(m1, 1.0f);
     125        LOLUNIT_ASSERT_EQUAL(m2, 2.0f);
     126        LOLUNIT_ASSERT_EQUAL(m3, 0.5f);
     127        LOLUNIT_ASSERT_EQUAL(m4, 1.0f);
     128    }
    113129};
    114130
Note: See TracChangeset for help on using the changeset viewer.