Changeset 1017


Ignore:
Timestamp:
Oct 12, 2011, 12:01:44 AM (9 years ago)
Author:
sam
Message:

core: fix an accuracy error in real::re() and real::sqrt() introduced in
the 16-to-32-bit refactoring.

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/real.cpp

    r1016 r1017  
    2323{
    2424
    25 static int const BIGIT_BITS = 32;
    26 
    2725real::real(float f) { *this = (double)f; }
    2826real::real(int i) { *this = (double)i; }
     
    501499    ret.m_signexp |= (exponent - 1) & 0x7fffffffu;
    502500
    503     /* FIXME: log2(BIGITS) steps of Newton-Raphson seems to be enough for
     501    /* FIXME: 1+log2(BIGITS) steps of Newton-Raphson seems to be enough for
    504502     * convergence, but this hasn't been checked seriously. */
    505     for (int i = 2; i < real::BIGITS; i *= 2)
     503    for (int i = 1; i <= real::BIGITS; i *= 2)
    506504        ret = ret * (real::R_2 - ret * x);
    507505
     
    545543    ret.m_signexp |= exponent & 0x7fffffffu;
    546544
    547     /* FIXME: log2(BIGITS) steps of Newton-Raphson seems to be enough for
     545    /* FIXME: 1+log2(BIGITS) steps of Newton-Raphson seems to be enough for
    548546     * convergence, but this hasn't been checked seriously. */
    549     for (int i = 2; i < real::BIGITS; i *= 2)
     547    for (int i = 1; i <= real::BIGITS; i *= 2)
    550548    {
    551549        ret = ret * (real::R_3 - ret * ret * x);
     
    669667        if (exponent <= 0)
    670668            ret.m_mantissa[i] = 0;
    671         else if (exponent < BIGIT_BITS)
    672             ret.m_mantissa[i] &= ~((1 << (BIGIT_BITS - exponent)) - 1);
    673 
    674         exponent -= BIGIT_BITS;
     669        else if (exponent < real::BIGIT_BITS)
     670            ret.m_mantissa[i] &= ~((1 << (real::BIGIT_BITS - exponent)) - 1);
     671
     672        exponent -= real::BIGIT_BITS;
    675673    }
    676674
  • trunk/src/real.h

    r1016 r1017  
    105105    static real const R_SQRT1_2;
    106106
    107 private:
    108107    /* XXX: changing this requires tuning real::fres (the number of
    109108     * Newton-Raphson iterations) and real::print (the number of printed
    110109     * digits) */
    111110    static int const BIGITS = 16;
     111    static int const BIGIT_BITS = 32;
    112112
     113private:
    113114    uint32_t m_size;
    114115    uint32_t m_signexp;
  • trunk/test/unit/real.cpp

    r1003 r1017  
    181181    }
    182182
    183     LOLUNIT_TEST(Division)
    184     {
    185         real a1(1.0f);
    186         real a2(2.0f);
    187 
    188         float m1 = a1 / a1;
    189         float m2 = a2 / a1;
    190         float m3 = a1 / a2;
    191         float m4 = a2 / a2;
    192         float m5 = a1 / -a2;
     183    LOLUNIT_TEST(ExactDivision)
     184    {
     185        float m1 = real::R_1 / real::R_1;
     186        float m2 = real::R_2 / real::R_1;
     187        float m3 = real::R_1 / real::R_2;
     188        float m4 = real::R_2 / real::R_2;
     189        float m5 = real::R_1 / -real::R_2;
    193190
    194191        LOLUNIT_ASSERT_EQUAL(m1, 1.0f);
     
    197194        LOLUNIT_ASSERT_EQUAL(m4, 1.0f);
    198195        LOLUNIT_ASSERT_EQUAL(m5, -0.5f);
     196    }
     197
     198    LOLUNIT_TEST(InexactDivision)
     199    {
     200        /* 1 / 3 * 3 should be close to 1... check that it does not differ
     201         * by more than 2^-k where k is the number of bits in the mantissa. */
     202        real a = real::R_1 / real::R_3 * real::R_3;
     203        real b = (real::R_1 - a) << (real::BIGITS * real::BIGIT_BITS);
     204
     205        LOLUNIT_ASSERT_LEQUAL((double)fabs(b), 1.0);
    199206    }
    200207
Note: See TracChangeset for help on using the changeset viewer.