Changeset 1122


Ignore:
Timestamp:
Jan 2, 2012, 4:14:23 AM (6 years ago)
Author:
sam
Message:

real: get rid of <<= and >>= operators; we can use ldexp() instead. As a
bonus, multiplication or division by a power of two will be optimised as
a shift of the exponent.

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/lol/math/real.h

    r1117 r1122  
    5656    real const &operator /=(real const &x);
    5757
    58     real operator <<(int x) const;
    59     real operator >>(int x) const;
    60     real const &operator <<=(int x);
    61     real const &operator >>=(int x);
    62 
    6358    bool operator ==(real const &x) const;
    6459    bool operator !=(real const &x) const;
     
    113108    void print(int ndigits = 150) const;
    114109
     110    /* Additional operators using base C++ types */
     111#define __LOL_REAL_OP_HELPER_GENERIC(op, type) \
     112    inline real operator op(type x) const { return *this op (real)x; } \
     113    inline real const &operator op##=(type x) { return *this = (*this op x); }
     114#define __LOL_REAL_OP_HELPER_FASTMULDIV(op, type) \
     115    inline real operator op(type x) const \
     116    { \
     117        real tmp = *this; return tmp op##= x; \
     118    } \
     119    inline real const &operator op##=(type x) \
     120    { \
     121        /* If multiplying or dividing by a power of two, take a shortcut */ \
     122        if ((m_signexp << 1) && x && !(x & (x - 1))) \
     123        { \
     124            while (x >>= 1) \
     125                m_signexp += 1 op 2 - 1; /* 1 if op is *, -1 if op is / */ \
     126        } \
     127        else \
     128            *this = *this op (real)x; \
     129        return *this; \
     130    }
     131#define __LOL_REAL_OP_HELPER_INT(type) \
     132    __LOL_REAL_OP_HELPER_GENERIC(+, type) \
     133    __LOL_REAL_OP_HELPER_GENERIC(-, type) \
     134    __LOL_REAL_OP_HELPER_FASTMULDIV(*, type) \
     135    __LOL_REAL_OP_HELPER_FASTMULDIV(/, type)
     136#define __LOL_REAL_OP_HELPER_FLOAT(type) \
     137    __LOL_REAL_OP_HELPER_GENERIC(+, type) \
     138    __LOL_REAL_OP_HELPER_GENERIC(-, type) \
     139    __LOL_REAL_OP_HELPER_GENERIC(*, type) \
     140    __LOL_REAL_OP_HELPER_GENERIC(/, type)
     141
     142    __LOL_REAL_OP_HELPER_INT(int)
     143    __LOL_REAL_OP_HELPER_INT(unsigned int)
     144    __LOL_REAL_OP_HELPER_FLOAT(float)
     145    __LOL_REAL_OP_HELPER_FLOAT(double)
     146
     147    /* Constants */
    115148    static real const R_0;
    116149    static real const R_1;
  • trunk/src/lol/math/remez.h

    r1119 r1122  
    3737        m_func = func;
    3838        m_weight = weight;
    39         m_k1 = (b + a) >> 1;
    40         m_k2 = (b - a) >> 1;
     39        m_k1 = (b + a) / 2;
     40        m_k2 = (b - a) / 2;
    4141        m_invk2 = re(m_k2);
    4242        m_invk1 = -m_k1 * m_invk2;
     
    140140            right.error = ChebyEval(right.value) - Value(right.value);
    141141
    142             static T limit = (T)1 >> 500;
     142            static T limit = ldexp((T)1, -500);
    143143            static T zero = (T)0;
    144144            while (fabs(left.value - right.value) > limit)
    145145            {
    146                 mid.value = (left.value + right.value) >> 1;
     146                mid.value = (left.value + right.value) / 2;
    147147                mid.error = ChebyEval(mid.value) - Value(mid.value);
    148148
     
    177177            for (;;)
    178178            {
    179                 T c = a, delta = (b - a) >> 3;
     179                T c = a, delta = (b - a) / 8;
    180180                T maxerror = 0;
    181181                T maxweight = 0;
     
    202202                    if (e > final)
    203203                        final = e;
    204                     control[i] = (a + b) >> 1;
     204                    control[i] = (a + b) / 2;
    205205                    break;
    206206                }
  • trunk/src/real.cpp

    r1121 r1122  
    464464}
    465465
    466 real real::operator <<(int x) const
    467 {
    468     real tmp = *this;
    469     return tmp <<= x;
    470 }
    471 
    472 real real::operator >>(int x) const
    473 {
    474     real tmp = *this;
    475     return tmp >>= x;
    476 }
    477 
    478 real const &real::operator <<=(int x)
    479 {
    480     if (m_signexp << 1)
    481         m_signexp += x;
    482     return *this;
    483 }
    484 
    485 real const &real::operator >>=(int x)
    486 {
    487     if (m_signexp << 1)
    488         m_signexp -= x;
    489     return *this;
    490 }
    491 
    492466bool real::operator ==(real const &x) const
    493467{
     
    680654    {
    681655        static real third = re(real::R_3);
    682         ret = third * (x / (ret * ret) + (ret << 1));
     656        ret = third * (x / (ret * ret) + (ret / 2));
    683657    }
    684658
     
    697671    {
    698672        /* Odd integer exponent */
    699         if (y == (round(y >> 1) << 1))
     673        if (y == (round(y / 2) * 2))
    700674            return exp(y * log(-x));
    701675
     
    718692    int a = ceilf(logf(2) / logf(2 * M_PI) * real::BIGITS * real::BIGIT_BITS);
    719693
    720     real ret = sqrt(real::R_PI << 1);
     694    real ret = sqrt(real::R_PI * 2);
    721695    real fact_k_1 = real::R_1;
    722696
     
    730704    }
    731705
    732     ret *= pow(x + (real)(a - 1), x - (real::R_1 >> 1));
     706    ret *= pow(x + (real)(a - 1), x - (real::R_1 / 2));
    733707    ret *= exp(-x - (real)(a - 1));
    734708
     
    773747    }
    774748
    775     return z * (sum << 2);
     749    return z * sum * 4;
    776750}
    777751
     
    876850     * accuracy near zero. We only use this identity for |x|>0.5. If
    877851     * |x|<=0.5, we compute exp(x)-1 and exp(-x)-1 instead. */
    878     bool near_zero = (fabs(x) < real::R_1 >> 1);
     852    bool near_zero = (fabs(x) < real::R_1 / 2);
    879853    real x1 = near_zero ? fast_exp_sub(x, real::R_1) : exp(x);
    880854    real x2 = near_zero ? fast_exp_sub(-x, real::R_1) : exp(-x);
    881     return (x1 - x2) >> 1;
     855    return (x1 - x2) / 2;
    882856}
    883857
     
    885859{
    886860    /* See sinh() for the strategy here */
    887     bool near_zero = (fabs(x) < real::R_1 >> 1);
     861    bool near_zero = (fabs(x) < real::R_1 / 2);
    888862    real x1 = near_zero ? fast_exp_sub(x, real::R_1) : exp(x);
    889863    real x2 = near_zero ? fast_exp_sub(-x, real::R_1) : exp(-x);
     
    896870    /* No need to worry about accuracy here; maybe the last bit is slightly
    897871     * off, but that's about it. */
    898     return (exp(x) + exp(-x)) >> 1;
     872    return (exp(x) + exp(-x)) / 2;
    899873}
    900874
     
    990964        return -round(-x);
    991965
    992     return floor(x + (real::R_1 >> 1));
     966    return floor(x + (real::R_1 / 2));
    993967}
    994968
     
    1009983    bool switch_sign = x.m_signexp & 0x80000000u;
    1010984
    1011     real absx = fmod(fabs(x), real::R_PI << 1);
     985    real absx = fmod(fabs(x), real::R_PI * 2);
    1012986    if (absx > real::R_PI)
    1013987    {
     
    10761050     * lose the precision around x=1. */
    10771051    real absx = fabs(x);
    1078     bool around_zero = (absx < (real::R_1 >> 1));
     1052    bool around_zero = (absx < (real::R_1 / 2));
    10791053
    10801054    if (!around_zero)
    1081         absx = sqrt((real::R_1 - absx) >> 1);
     1055        absx = sqrt((real::R_1 - absx) / 2);
    10821056
    10831057    real ret = absx, xn = absx, x2 = absx * absx, fact1 = 2, fact2 = 1;
     
    10861060        xn *= x2;
    10871061        real mul = (real)(2 * i + 1);
    1088         real newret = ret + ((fact1 * xn / (mul * fact2)) >> (i * 2));
     1062        real newret = ret + ldexp(fact1 * xn / (mul * fact2), -2 * i);
    10891063        if (newret == ret)
    10901064            break;
     
    11031077        real adjust = is_negative ? real::R_PI : real::R_0;
    11041078        if (is_asin)
    1105             ret = real::R_PI_2 - adjust - (ret << 1);
     1079            ret = real::R_PI_2 - adjust - ret * 2;
    11061080        else
    1107             ret = adjust + (ret << 1);
     1081            ret = adjust + ret * 2;
    11081082    }
    11091083
     
    11471121    real absx = fabs(x);
    11481122
    1149     if (absx < (real::R_1 >> 1))
     1123    if (absx < (real::R_1 / 2))
    11501124    {
    11511125        real ret = x, xn = x, mx2 = -x * x;
     
    11631137    real ret = 0;
    11641138
    1165     if (absx < (real::R_3 >> 1))
     1139    if (absx < (real::R_3 / 2))
    11661140    {
    11671141        real y = real::R_1 - absx;
     
    11691143        for (int i = 0; ; i += 2)
    11701144        {
    1171             real newret = ret + ((yn / (real)(2 * i + 1)) >> (i + 1));
     1145            real newret = ret + ldexp(yn / (real)(2 * i + 1), -i - 1);
    11721146            yn *= y;
    1173             newret += (yn / (real)(2 * i + 2)) >> (i + 1);
     1147            newret += ldexp(yn / (real)(2 * i + 2), -i - 1);
    11741148            yn *= y;
    1175             newret += (yn / (real)(2 * i + 3)) >> (i + 2);
     1149            newret += ldexp(yn / (real)(2 * i + 3), -i - 2);
    11761150            if (newret == ret)
    11771151                break;
     
    11831157    else if (absx < real::R_2)
    11841158    {
    1185         real y = (absx - real::R_SQRT3) >> 1;
     1159        real y = (absx - real::R_SQRT3) / 2;
    11861160        real yn = y, my2 = -y * y;
    11871161        for (int i = 1; ; i += 6)
    11881162        {
    1189             real newret = ret + ((yn / (real)i) >> 1);
     1163            real newret = ret + ((yn / (real)i) / 2);
    11901164            yn *= y;
    1191             newret -= (real::R_SQRT3 >> 1) * yn / (real)(i + 1);
     1165            newret -= (real::R_SQRT3 / 2) * yn / (real)(i + 1);
    11921166            yn *= y;
    11931167            newret += yn / (real)(i + 2);
    11941168            yn *= y;
    1195             newret -= (real::R_SQRT3 >> 1) * yn / (real)(i + 3);
     1169            newret -= (real::R_SQRT3 / 2) * yn / (real)(i + 3);
    11961170            yn *= y;
    1197             newret += (yn / (real)(i + 4)) >> 1;
     1171            newret += (yn / (real)(i + 4)) / 2;
    11981172            if (newret == ret)
    11991173                break;
     
    13301304real const real::R_E        = exp(R_1);
    13311305real const real::R_PI       = fast_pi();
    1332 real const real::R_PI_2     = R_PI >> 1;
     1306real const real::R_PI_2     = R_PI / 2;
    13331307real const real::R_PI_3     = R_PI / R_3;
    1334 real const real::R_PI_4     = R_PI >> 2;
     1308real const real::R_PI_4     = R_PI / 4;
    13351309real const real::R_1_PI     = re(R_PI);
    1336 real const real::R_2_PI     = R_1_PI << 1;
    1337 real const real::R_2_SQRTPI = re(sqrt(R_PI)) << 1;
     1310real const real::R_2_PI     = R_1_PI * 2;
     1311real const real::R_2_SQRTPI = re(sqrt(R_PI)) * 2;
    13381312real const real::R_SQRT2    = sqrt(R_2);
    13391313real const real::R_SQRT3    = sqrt(R_3);
    1340 real const real::R_SQRT1_2  = R_SQRT2 >> 1;
     1314real const real::R_SQRT1_2  = R_SQRT2 / 2;
    13411315
    13421316} /* namespace lol */
  • trunk/test/unit/real.cpp

    r1116 r1122  
    222222         * by more than 2^-k where k is the number of bits in the mantissa. */
    223223        real a = real::R_1 / real::R_3 * real::R_3;
    224         real b = (real::R_1 - a) << (real::BIGITS * real::BIGIT_BITS);
     224        real b = ldexp(real::R_1 - a, real::BIGITS * real::BIGIT_BITS);
    225225
    226226        LOLUNIT_ASSERT_LEQUAL((double)fabs(b), 1.0);
    227227    }
    228228
    229     LOLUNIT_TEST(Shift)
     229    LOLUNIT_TEST(LoadExp)
    230230    {
    231231        real a1(1.5);
     
    233233        real a3(0.0);
    234234
    235         LOLUNIT_ASSERT_EQUAL((double)(a1 >> 7), 0.01171875);
    236         LOLUNIT_ASSERT_EQUAL((double)(a1 >> -7), 192.0);
    237         LOLUNIT_ASSERT_EQUAL((double)(a1 << 7), 192.0);
    238         LOLUNIT_ASSERT_EQUAL((double)(a1 << -7), 0.01171875);
    239 
    240         LOLUNIT_ASSERT_EQUAL((double)(a2 >> 7), -0.01171875);
    241         LOLUNIT_ASSERT_EQUAL((double)(a2 >> -7), -192.0);
    242         LOLUNIT_ASSERT_EQUAL((double)(a2 << 7), -192.0);
    243         LOLUNIT_ASSERT_EQUAL((double)(a2 << -7), -0.01171875);
    244 
    245         LOLUNIT_ASSERT_EQUAL((double)(a3 >> 7), 0.0);
    246         LOLUNIT_ASSERT_EQUAL((double)(a3 >> -7), 0.0);
    247         LOLUNIT_ASSERT_EQUAL((double)(a3 << 7), 0.0);
    248         LOLUNIT_ASSERT_EQUAL((double)(a3 << -7), 0.0);
     235        LOLUNIT_ASSERT_EQUAL((double)ldexp(a1, 7), 192.0);
     236        LOLUNIT_ASSERT_EQUAL((double)ldexp(a1, -7), 0.01171875);
     237
     238        LOLUNIT_ASSERT_EQUAL((double)ldexp(a2, 7), -192.0);
     239        LOLUNIT_ASSERT_EQUAL((double)ldexp(a2, -7), -0.01171875);
     240
     241        LOLUNIT_ASSERT_EQUAL((double)ldexp(a3, 7), 0.0);
     242        LOLUNIT_ASSERT_EQUAL((double)ldexp(a3, -7), 0.0);
    249243    }
    250244
Note: See TracChangeset for help on using the changeset viewer.