# Changeset 1122Tweet

Ignore:
Timestamp:
Jan 2, 2012, 4:14:23 AM (6 years ago)
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

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

 r1117 real const &operator /=(real const &x); real operator <<(int x) const; real operator >>(int x) const; real const &operator <<=(int x); real const &operator >>=(int x); bool operator ==(real const &x) const; bool operator !=(real const &x) const; void print(int ndigits = 150) const; /* Additional operators using base C++ types */ #define __LOL_REAL_OP_HELPER_GENERIC(op, type) \ inline real operator op(type x) const { return *this op (real)x; } \ inline real const &operator op##=(type x) { return *this = (*this op x); } #define __LOL_REAL_OP_HELPER_FASTMULDIV(op, type) \ inline real operator op(type x) const \ { \ real tmp = *this; return tmp op##= x; \ } \ inline real const &operator op##=(type x) \ { \ /* If multiplying or dividing by a power of two, take a shortcut */ \ if ((m_signexp << 1) && x && !(x & (x - 1))) \ { \ while (x >>= 1) \ m_signexp += 1 op 2 - 1; /* 1 if op is *, -1 if op is / */ \ } \ else \ *this = *this op (real)x; \ return *this; \ } #define __LOL_REAL_OP_HELPER_INT(type) \ __LOL_REAL_OP_HELPER_GENERIC(+, type) \ __LOL_REAL_OP_HELPER_GENERIC(-, type) \ __LOL_REAL_OP_HELPER_FASTMULDIV(*, type) \ __LOL_REAL_OP_HELPER_FASTMULDIV(/, type) #define __LOL_REAL_OP_HELPER_FLOAT(type) \ __LOL_REAL_OP_HELPER_GENERIC(+, type) \ __LOL_REAL_OP_HELPER_GENERIC(-, type) \ __LOL_REAL_OP_HELPER_GENERIC(*, type) \ __LOL_REAL_OP_HELPER_GENERIC(/, type) __LOL_REAL_OP_HELPER_INT(int) __LOL_REAL_OP_HELPER_INT(unsigned int) __LOL_REAL_OP_HELPER_FLOAT(float) __LOL_REAL_OP_HELPER_FLOAT(double) /* Constants */ static real const R_0; static real const R_1;
• ## trunk/src/lol/math/remez.h

 r1119 m_func = func; m_weight = weight; m_k1 = (b + a) >> 1; m_k2 = (b - a) >> 1; m_k1 = (b + a) / 2; m_k2 = (b - a) / 2; m_invk2 = re(m_k2); m_invk1 = -m_k1 * m_invk2; right.error = ChebyEval(right.value) - Value(right.value); static T limit = (T)1 >> 500; static T limit = ldexp((T)1, -500); static T zero = (T)0; while (fabs(left.value - right.value) > limit) { mid.value = (left.value + right.value) >> 1; mid.value = (left.value + right.value) / 2; mid.error = ChebyEval(mid.value) - Value(mid.value); for (;;) { T c = a, delta = (b - a) >> 3; T c = a, delta = (b - a) / 8; T maxerror = 0; T maxweight = 0; if (e > final) final = e; control[i] = (a + b) >> 1; control[i] = (a + b) / 2; break; }
• ## trunk/src/real.cpp

 r1121 } real real::operator <<(int x) const { real tmp = *this; return tmp <<= x; } real real::operator >>(int x) const { real tmp = *this; return tmp >>= x; } real const &real::operator <<=(int x) { if (m_signexp << 1) m_signexp += x; return *this; } real const &real::operator >>=(int x) { if (m_signexp << 1) m_signexp -= x; return *this; } bool real::operator ==(real const &x) const { { static real third = re(real::R_3); ret = third * (x / (ret * ret) + (ret << 1)); ret = third * (x / (ret * ret) + (ret / 2)); } { /* Odd integer exponent */ if (y == (round(y >> 1) << 1)) if (y == (round(y / 2) * 2)) return exp(y * log(-x)); int a = ceilf(logf(2) / logf(2 * M_PI) * real::BIGITS * real::BIGIT_BITS); real ret = sqrt(real::R_PI << 1); real ret = sqrt(real::R_PI * 2); real fact_k_1 = real::R_1; } ret *= pow(x + (real)(a - 1), x - (real::R_1 >> 1)); ret *= pow(x + (real)(a - 1), x - (real::R_1 / 2)); ret *= exp(-x - (real)(a - 1)); } return z * (sum << 2); return z * sum * 4; } * accuracy near zero. We only use this identity for |x|>0.5. If * |x|<=0.5, we compute exp(x)-1 and exp(-x)-1 instead. */ bool near_zero = (fabs(x) < real::R_1 >> 1); bool near_zero = (fabs(x) < real::R_1 / 2); real x1 = near_zero ? fast_exp_sub(x, real::R_1) : exp(x); real x2 = near_zero ? fast_exp_sub(-x, real::R_1) : exp(-x); return (x1 - x2) >> 1; return (x1 - x2) / 2; } { /* See sinh() for the strategy here */ bool near_zero = (fabs(x) < real::R_1 >> 1); bool near_zero = (fabs(x) < real::R_1 / 2); real x1 = near_zero ? fast_exp_sub(x, real::R_1) : exp(x); real x2 = near_zero ? fast_exp_sub(-x, real::R_1) : exp(-x); /* No need to worry about accuracy here; maybe the last bit is slightly * off, but that's about it. */ return (exp(x) + exp(-x)) >> 1; return (exp(x) + exp(-x)) / 2; } return -round(-x); return floor(x + (real::R_1 >> 1)); return floor(x + (real::R_1 / 2)); } bool switch_sign = x.m_signexp & 0x80000000u; real absx = fmod(fabs(x), real::R_PI << 1); real absx = fmod(fabs(x), real::R_PI * 2); if (absx > real::R_PI) { * lose the precision around x=1. */ real absx = fabs(x); bool around_zero = (absx < (real::R_1 >> 1)); bool around_zero = (absx < (real::R_1 / 2)); if (!around_zero) absx = sqrt((real::R_1 - absx) >> 1); absx = sqrt((real::R_1 - absx) / 2); real ret = absx, xn = absx, x2 = absx * absx, fact1 = 2, fact2 = 1; xn *= x2; real mul = (real)(2 * i + 1); real newret = ret + ((fact1 * xn / (mul * fact2)) >> (i * 2)); real newret = ret + ldexp(fact1 * xn / (mul * fact2), -2 * i); if (newret == ret) break; real adjust = is_negative ? real::R_PI : real::R_0; if (is_asin) ret = real::R_PI_2 - adjust - (ret << 1); ret = real::R_PI_2 - adjust - ret * 2; else ret = adjust + (ret << 1); ret = adjust + ret * 2; } real absx = fabs(x); if (absx < (real::R_1 >> 1)) if (absx < (real::R_1 / 2)) { real ret = x, xn = x, mx2 = -x * x; real ret = 0; if (absx < (real::R_3 >> 1)) if (absx < (real::R_3 / 2)) { real y = real::R_1 - absx; for (int i = 0; ; i += 2) { real newret = ret + ((yn / (real)(2 * i + 1)) >> (i + 1)); real newret = ret + ldexp(yn / (real)(2 * i + 1), -i - 1); yn *= y; newret += (yn / (real)(2 * i + 2)) >> (i + 1); newret += ldexp(yn / (real)(2 * i + 2), -i - 1); yn *= y; newret += (yn / (real)(2 * i + 3)) >> (i + 2); newret += ldexp(yn / (real)(2 * i + 3), -i - 2); if (newret == ret) break; else if (absx < real::R_2) { real y = (absx - real::R_SQRT3) >> 1; real y = (absx - real::R_SQRT3) / 2; real yn = y, my2 = -y * y; for (int i = 1; ; i += 6) { real newret = ret + ((yn / (real)i) >> 1); real newret = ret + ((yn / (real)i) / 2); yn *= y; newret -= (real::R_SQRT3 >> 1) * yn / (real)(i + 1); newret -= (real::R_SQRT3 / 2) * yn / (real)(i + 1); yn *= y; newret += yn / (real)(i + 2); yn *= y; newret -= (real::R_SQRT3 >> 1) * yn / (real)(i + 3); newret -= (real::R_SQRT3 / 2) * yn / (real)(i + 3); yn *= y; newret += (yn / (real)(i + 4)) >> 1; newret += (yn / (real)(i + 4)) / 2; if (newret == ret) break; real const real::R_E        = exp(R_1); real const real::R_PI       = fast_pi(); real const real::R_PI_2     = R_PI >> 1; real const real::R_PI_2     = R_PI / 2; real const real::R_PI_3     = R_PI / R_3; real const real::R_PI_4     = R_PI >> 2; real const real::R_PI_4     = R_PI / 4; real const real::R_1_PI     = re(R_PI); real const real::R_2_PI     = R_1_PI << 1; real const real::R_2_SQRTPI = re(sqrt(R_PI)) << 1; real const real::R_2_PI     = R_1_PI * 2; real const real::R_2_SQRTPI = re(sqrt(R_PI)) * 2; real const real::R_SQRT2    = sqrt(R_2); real const real::R_SQRT3    = sqrt(R_3); real const real::R_SQRT1_2  = R_SQRT2 >> 1; real const real::R_SQRT1_2  = R_SQRT2 / 2; } /* namespace lol */
• ## trunk/test/unit/real.cpp

 r1116 * by more than 2^-k where k is the number of bits in the mantissa. */ real a = real::R_1 / real::R_3 * real::R_3; real b = (real::R_1 - a) << (real::BIGITS * real::BIGIT_BITS); real b = ldexp(real::R_1 - a, real::BIGITS * real::BIGIT_BITS); LOLUNIT_ASSERT_LEQUAL((double)fabs(b), 1.0); } LOLUNIT_TEST(Shift) LOLUNIT_TEST(LoadExp) { real a1(1.5); real a3(0.0); LOLUNIT_ASSERT_EQUAL((double)(a1 >> 7), 0.01171875); LOLUNIT_ASSERT_EQUAL((double)(a1 >> -7), 192.0); LOLUNIT_ASSERT_EQUAL((double)(a1 << 7), 192.0); LOLUNIT_ASSERT_EQUAL((double)(a1 << -7), 0.01171875); LOLUNIT_ASSERT_EQUAL((double)(a2 >> 7), -0.01171875); LOLUNIT_ASSERT_EQUAL((double)(a2 >> -7), -192.0); LOLUNIT_ASSERT_EQUAL((double)(a2 << 7), -192.0); LOLUNIT_ASSERT_EQUAL((double)(a2 << -7), -0.01171875); LOLUNIT_ASSERT_EQUAL((double)(a3 >> 7), 0.0); LOLUNIT_ASSERT_EQUAL((double)(a3 >> -7), 0.0); LOLUNIT_ASSERT_EQUAL((double)(a3 << 7), 0.0); LOLUNIT_ASSERT_EQUAL((double)(a3 << -7), 0.0); LOLUNIT_ASSERT_EQUAL((double)ldexp(a1, 7), 192.0); LOLUNIT_ASSERT_EQUAL((double)ldexp(a1, -7), 0.01171875); LOLUNIT_ASSERT_EQUAL((double)ldexp(a2, 7), -192.0); LOLUNIT_ASSERT_EQUAL((double)ldexp(a2, -7), -0.01171875); LOLUNIT_ASSERT_EQUAL((double)ldexp(a3, 7), 0.0); LOLUNIT_ASSERT_EQUAL((double)ldexp(a3, -7), 0.0); }
Note: See TracChangeset for help on using the changeset viewer.