Changeset 1122
 Timestamp:
 Jan 2, 2012, 4:14:23 AM (7 years ago)
 Location:
 trunk
 Files:

 4 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/lol/math/real.h
r1117 r1122 56 56 real const &operator /=(real const &x); 57 57 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 63 58 bool operator ==(real const &x) const; 64 59 bool operator !=(real const &x) const; … … 113 108 void print(int ndigits = 150) const; 114 109 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 */ 115 148 static real const R_0; 116 149 static real const R_1; 
trunk/src/lol/math/remez.h
r1119 r1122 37 37 m_func = func; 38 38 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; 41 41 m_invk2 = re(m_k2); 42 42 m_invk1 = m_k1 * m_invk2; … … 140 140 right.error = ChebyEval(right.value)  Value(right.value); 141 141 142 static T limit = (T)1 >> 500;142 static T limit = ldexp((T)1, 500); 143 143 static T zero = (T)0; 144 144 while (fabs(left.value  right.value) > limit) 145 145 { 146 mid.value = (left.value + right.value) >> 1;146 mid.value = (left.value + right.value) / 2; 147 147 mid.error = ChebyEval(mid.value)  Value(mid.value); 148 148 … … 177 177 for (;;) 178 178 { 179 T c = a, delta = (b  a) >> 3;179 T c = a, delta = (b  a) / 8; 180 180 T maxerror = 0; 181 181 T maxweight = 0; … … 202 202 if (e > final) 203 203 final = e; 204 control[i] = (a + b) >> 1;204 control[i] = (a + b) / 2; 205 205 break; 206 206 } 
trunk/src/real.cpp
r1121 r1122 464 464 } 465 465 466 real real::operator <<(int x) const467 {468 real tmp = *this;469 return tmp <<= x;470 }471 472 real real::operator >>(int x) const473 {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 492 466 bool real::operator ==(real const &x) const 493 467 { … … 680 654 { 681 655 static real third = re(real::R_3); 682 ret = third * (x / (ret * ret) + (ret << 1));656 ret = third * (x / (ret * ret) + (ret / 2)); 683 657 } 684 658 … … 697 671 { 698 672 /* Odd integer exponent */ 699 if (y == (round(y >> 1) << 1))673 if (y == (round(y / 2) * 2)) 700 674 return exp(y * log(x)); 701 675 … … 718 692 int a = ceilf(logf(2) / logf(2 * M_PI) * real::BIGITS * real::BIGIT_BITS); 719 693 720 real ret = sqrt(real::R_PI << 1);694 real ret = sqrt(real::R_PI * 2); 721 695 real fact_k_1 = real::R_1; 722 696 … … 730 704 } 731 705 732 ret *= pow(x + (real)(a  1), x  (real::R_1 >> 1));706 ret *= pow(x + (real)(a  1), x  (real::R_1 / 2)); 733 707 ret *= exp(x  (real)(a  1)); 734 708 … … 773 747 } 774 748 775 return z * (sum << 2);749 return z * sum * 4; 776 750 } 777 751 … … 876 850 * accuracy near zero. We only use this identity for x>0.5. If 877 851 * 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); 879 853 real x1 = near_zero ? fast_exp_sub(x, real::R_1) : exp(x); 880 854 real x2 = near_zero ? fast_exp_sub(x, real::R_1) : exp(x); 881 return (x1  x2) >> 1;855 return (x1  x2) / 2; 882 856 } 883 857 … … 885 859 { 886 860 /* 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); 888 862 real x1 = near_zero ? fast_exp_sub(x, real::R_1) : exp(x); 889 863 real x2 = near_zero ? fast_exp_sub(x, real::R_1) : exp(x); … … 896 870 /* No need to worry about accuracy here; maybe the last bit is slightly 897 871 * off, but that's about it. */ 898 return (exp(x) + exp(x)) >> 1;872 return (exp(x) + exp(x)) / 2; 899 873 } 900 874 … … 990 964 return round(x); 991 965 992 return floor(x + (real::R_1 >> 1));966 return floor(x + (real::R_1 / 2)); 993 967 } 994 968 … … 1009 983 bool switch_sign = x.m_signexp & 0x80000000u; 1010 984 1011 real absx = fmod(fabs(x), real::R_PI << 1);985 real absx = fmod(fabs(x), real::R_PI * 2); 1012 986 if (absx > real::R_PI) 1013 987 { … … 1076 1050 * lose the precision around x=1. */ 1077 1051 real absx = fabs(x); 1078 bool around_zero = (absx < (real::R_1 >> 1));1052 bool around_zero = (absx < (real::R_1 / 2)); 1079 1053 1080 1054 if (!around_zero) 1081 absx = sqrt((real::R_1  absx) >> 1);1055 absx = sqrt((real::R_1  absx) / 2); 1082 1056 1083 1057 real ret = absx, xn = absx, x2 = absx * absx, fact1 = 2, fact2 = 1; … … 1086 1060 xn *= x2; 1087 1061 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); 1089 1063 if (newret == ret) 1090 1064 break; … … 1103 1077 real adjust = is_negative ? real::R_PI : real::R_0; 1104 1078 if (is_asin) 1105 ret = real::R_PI_2  adjust  (ret << 1);1079 ret = real::R_PI_2  adjust  ret * 2; 1106 1080 else 1107 ret = adjust + (ret << 1);1081 ret = adjust + ret * 2; 1108 1082 } 1109 1083 … … 1147 1121 real absx = fabs(x); 1148 1122 1149 if (absx < (real::R_1 >> 1))1123 if (absx < (real::R_1 / 2)) 1150 1124 { 1151 1125 real ret = x, xn = x, mx2 = x * x; … … 1163 1137 real ret = 0; 1164 1138 1165 if (absx < (real::R_3 >> 1))1139 if (absx < (real::R_3 / 2)) 1166 1140 { 1167 1141 real y = real::R_1  absx; … … 1169 1143 for (int i = 0; ; i += 2) 1170 1144 { 1171 real newret = ret + ((yn / (real)(2 * i + 1)) >> (i + 1));1145 real newret = ret + ldexp(yn / (real)(2 * i + 1), i  1); 1172 1146 yn *= y; 1173 newret += (yn / (real)(2 * i + 2)) >> (i +1);1147 newret += ldexp(yn / (real)(2 * i + 2), i  1); 1174 1148 yn *= y; 1175 newret += (yn / (real)(2 * i + 3)) >> (i +2);1149 newret += ldexp(yn / (real)(2 * i + 3), i  2); 1176 1150 if (newret == ret) 1177 1151 break; … … 1183 1157 else if (absx < real::R_2) 1184 1158 { 1185 real y = (absx  real::R_SQRT3) >> 1;1159 real y = (absx  real::R_SQRT3) / 2; 1186 1160 real yn = y, my2 = y * y; 1187 1161 for (int i = 1; ; i += 6) 1188 1162 { 1189 real newret = ret + ((yn / (real)i) >> 1);1163 real newret = ret + ((yn / (real)i) / 2); 1190 1164 yn *= y; 1191 newret = (real::R_SQRT3 >> 1) * yn / (real)(i + 1);1165 newret = (real::R_SQRT3 / 2) * yn / (real)(i + 1); 1192 1166 yn *= y; 1193 1167 newret += yn / (real)(i + 2); 1194 1168 yn *= y; 1195 newret = (real::R_SQRT3 >> 1) * yn / (real)(i + 3);1169 newret = (real::R_SQRT3 / 2) * yn / (real)(i + 3); 1196 1170 yn *= y; 1197 newret += (yn / (real)(i + 4)) >> 1;1171 newret += (yn / (real)(i + 4)) / 2; 1198 1172 if (newret == ret) 1199 1173 break; … … 1330 1304 real const real::R_E = exp(R_1); 1331 1305 real const real::R_PI = fast_pi(); 1332 real const real::R_PI_2 = R_PI >> 1;1306 real const real::R_PI_2 = R_PI / 2; 1333 1307 real const real::R_PI_3 = R_PI / R_3; 1334 real const real::R_PI_4 = R_PI >> 2;1308 real const real::R_PI_4 = R_PI / 4; 1335 1309 real 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;1310 real const real::R_2_PI = R_1_PI * 2; 1311 real const real::R_2_SQRTPI = re(sqrt(R_PI)) * 2; 1338 1312 real const real::R_SQRT2 = sqrt(R_2); 1339 1313 real const real::R_SQRT3 = sqrt(R_3); 1340 real const real::R_SQRT1_2 = R_SQRT2 >> 1;1314 real const real::R_SQRT1_2 = R_SQRT2 / 2; 1341 1315 1342 1316 } /* namespace lol */ 
trunk/test/unit/real.cpp
r1116 r1122 222 222 * by more than 2^k where k is the number of bits in the mantissa. */ 223 223 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); 225 225 226 226 LOLUNIT_ASSERT_LEQUAL((double)fabs(b), 1.0); 227 227 } 228 228 229 LOLUNIT_TEST( Shift)229 LOLUNIT_TEST(LoadExp) 230 230 { 231 231 real a1(1.5); … … 233 233 real a3(0.0); 234 234 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); 249 243 } 250 244
Note: See TracChangeset
for help on using the changeset viewer.