Changeset 1893 for trunk/src/math/real.cpp
 Timestamp:
 Sep 9, 2012, 4:55:26 PM (8 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/math/real.cpp
r1667 r1893 37 37 namespace lol 38 38 { 39 40 /* 41 * First handle explicit specialisation of our templates. 42 * 43 * Initialisation order is not important because everything is 44 * done on demand, but here is the dependency list anyway: 45 *  fast_log() requires R_1 46 *  log() requires R_LN2 47 *  re() require R_2 48 *  exp() requires R_0, R_1, R_LN2 49 *  sqrt() requires R_3 50 */ 51 52 static real fast_log(real const &x); 53 static real fast_pi(); 54 55 #define LOL_CONSTANT_GETTER(name, value) \ 56 template<> real const& real::name() \ 57 { \ 58 static real const ret = value; \ 59 return ret; \ 60 } 61 62 LOL_CONSTANT_GETTER(R_0, (real)0.0); 63 LOL_CONSTANT_GETTER(R_1, (real)1.0); 64 LOL_CONSTANT_GETTER(R_2, (real)2.0); 65 LOL_CONSTANT_GETTER(R_3, (real)3.0); 66 LOL_CONSTANT_GETTER(R_10, (real)10.0); 67 68 LOL_CONSTANT_GETTER(R_LN2, fast_log(R_2())); 69 LOL_CONSTANT_GETTER(R_LN10, log(R_10())); 70 LOL_CONSTANT_GETTER(R_LOG2E, re(R_LN2())); 71 LOL_CONSTANT_GETTER(R_LOG10E, re(R_LN10())); 72 LOL_CONSTANT_GETTER(R_E, exp(R_1())); 73 LOL_CONSTANT_GETTER(R_PI, fast_pi()); 74 LOL_CONSTANT_GETTER(R_PI_2, R_PI() / 2); 75 LOL_CONSTANT_GETTER(R_PI_3, R_PI() / R_3()); 76 LOL_CONSTANT_GETTER(R_PI_4, R_PI() / 4); 77 LOL_CONSTANT_GETTER(R_1_PI, re(R_PI())); 78 LOL_CONSTANT_GETTER(R_2_PI, R_1_PI() * 2); 79 LOL_CONSTANT_GETTER(R_2_SQRTPI, re(sqrt(R_PI())) * 2); 80 LOL_CONSTANT_GETTER(R_SQRT2, sqrt(R_2())); 81 LOL_CONSTANT_GETTER(R_SQRT3, sqrt(R_3())); 82 LOL_CONSTANT_GETTER(R_SQRT1_2, R_SQRT2() / 2); 83 84 #undef LOL_CONSTANT_GETTER 85 86 /* 87 * Now carry on with the rest of the Real class. 88 */ 39 89 40 90 template<> real::Real() … … 190 240 191 241 if (exponent) 192 ret *= pow(R_10 , (real)exponent);242 ret *= pow(R_10(), (real)exponent); 193 243 194 244 if (negative) … … 602 652 * convergence, but this hasn't been checked seriously. */ 603 653 for (int i = 1; i <= real::BIGITS; i *= 2) 604 ret = ret * (real::R_2  ret * x);654 ret = ret * (real::R_2()  ret * x); 605 655 606 656 return ret; … … 647 697 for (int i = 1; i <= real::BIGITS; i *= 2) 648 698 { 649 ret = ret * (real::R_3  ret * ret * x);699 ret = ret * (real::R_3()  ret * ret * x); 650 700 ret.m_signexp; 651 701 } … … 684 734 for (int i = 1; i <= real::BIGITS; i *= 2) 685 735 { 686 static real third = re(real::R_3 );736 static real third = re(real::R_3()); 687 737 ret = third * (x / (ret * ret) + (ret / 2)); 688 738 } … … 694 744 { 695 745 if (!y) 696 return real::R_1 ;746 return real::R_1(); 697 747 if (!x) 698 return real::R_0 ;699 if (x > real::R_0 )748 return real::R_0(); 749 if (x > real::R_0()) 700 750 return exp(y * log(x)); 701 751 else /* x < 0 */ … … 710 760 711 761 /* FIXME: negative nth root */ 712 return real::R_0 ;762 return real::R_0(); 713 763 } 714 764 } … … 716 766 static real fast_fact(int x) 717 767 { 718 real ret = real::R_1 ;768 real ret = real::R_1(); 719 769 int i = 1, multiplier = 1, exponent = 0; 720 770 … … 749 799 int a = ceilf(logf(2) / logf(2 * M_PI) * real::BIGITS * real::BIGIT_BITS); 750 800 751 real ret = sqrt(real::R_PI * 2);752 real fact_k_1 = real::R_1 ;801 real ret = sqrt(real::R_PI() * 2); 802 real fact_k_1 = real::R_1(); 753 803 754 804 for (int k = 1; k < a; k++) … … 761 811 } 762 812 763 ret *= pow(x + (real)(a  1), x  (real::R_1 / 2));813 ret *= pow(x + (real)(a  1), x  (real::R_1() / 2)); 764 814 ret *= exp(x  (real)(a  1)); 765 815 … … 792 842 * sqrt() call. */ 793 843 real y = sqrt(x); 794 real z = (y  real::R_1 ) / (y + real::R_1), z2 = z * z, zn = z2;795 real sum = real::R_1 ;844 real z = (y  real::R_1()) / (y + real::R_1()), z2 = z * z, zn = z2; 845 real sum = real::R_1(); 796 846 797 847 for (int i = 3; ; i += 2) … … 819 869 } 820 870 tmp.m_signexp = (1 << 30)  1; 821 return (real)(int)(x.m_signexp  (1 << 30) + 1) * real::R_LN2 871 return (real)(int)(x.m_signexp  (1 << 30) + 1) * real::R_LN2() 822 872 + fast_log(tmp); 823 873 } … … 835 885 tmp.m_signexp = (1 << 30)  1; 836 886 return (real)(int)(x.m_signexp  (1 << 30) + 1) 837 + fast_log(tmp) * real::R_LOG2E ;887 + fast_log(tmp) * real::R_LOG2E(); 838 888 } 839 889 840 890 template<> real log10(real const &x) 841 891 { 842 return log(x) * real::R_LOG10E ;892 return log(x) * real::R_LOG10E(); 843 893 } 844 894 … … 849 899 * domain of validity. The argument y is used for cases where we 850 900 * don't want the leading 1 in the Taylor series. */ 851 real ret = real::R_1  y, xn = x;901 real ret = real::R_1()  y, xn = x; 852 902 int i = 1; 853 903 … … 882 932 * return x1 * 2^E0 883 933 */ 884 int e0 = x / real::R_LN2 ;885 real x0 = x  (real)e0 * real::R_LN2 ;886 real x1 = fast_exp_sub(x0, real::R_0 );934 int e0 = x / real::R_LN2(); 935 real x0 = x  (real)e0 * real::R_LN2(); 936 real x1 = fast_exp_sub(x0, real::R_0()); 887 937 x1.m_signexp += e0; 888 938 return x1; … … 894 944 int e0 = x; 895 945 real x0 = x  (real)e0; 896 real x1 = fast_exp_sub(x0 * real::R_LN2 , real::R_0);946 real x1 = fast_exp_sub(x0 * real::R_LN2(), real::R_0()); 897 947 x1.m_signexp += e0; 898 948 return x1; … … 904 954 * accuracy near zero. We only use this identity for x>0.5. If 905 955 * x<=0.5, we compute exp(x)1 and exp(x)1 instead. */ 906 bool near_zero = (fabs(x) < real::R_1 / 2);907 real x1 = near_zero ? fast_exp_sub(x, real::R_1 ) : exp(x);908 real x2 = near_zero ? fast_exp_sub(x, real::R_1 ) : exp(x);956 bool near_zero = (fabs(x) < real::R_1() / 2); 957 real x1 = near_zero ? fast_exp_sub(x, real::R_1()) : exp(x); 958 real x2 = near_zero ? fast_exp_sub(x, real::R_1()) : exp(x); 909 959 return (x1  x2) / 2; 910 960 } … … 913 963 { 914 964 /* See sinh() for the strategy here */ 915 bool near_zero = (fabs(x) < real::R_1 / 2);916 real x1 = near_zero ? fast_exp_sub(x, real::R_1 ) : exp(x);917 real x2 = near_zero ? fast_exp_sub(x, real::R_1 ) : exp(x);918 real x3 = near_zero ? x1 + x2 + real::R_2 : x1 + x2;965 bool near_zero = (fabs(x) < real::R_1() / 2); 966 real x1 = near_zero ? fast_exp_sub(x, real::R_1()) : exp(x); 967 real x2 = near_zero ? fast_exp_sub(x, real::R_1()) : exp(x); 968 real x3 = near_zero ? x1 + x2 + real::R_2() : x1 + x2; 919 969 return (x1  x2) / x3; 920 970 } … … 961 1011 template<> real ulp(real const &x) 962 1012 { 963 real ret = real::R_1 ;1013 real ret = real::R_1(); 964 1014 if (x) 965 1015 ret.m_signexp = x.m_signexp + 1  real::BIGITS * real::BIGIT_BITS; … … 995 1045 *  otherwise, if e is the exponent, clear all bits except the 996 1046 * first e. */ 997 if (x < real::R_0 )1047 if (x < real::R_0()) 998 1048 return ceil(x); 999 1049 if (!x) 1000 1050 return x; 1001 if (x < real::R_1 )1002 return real::R_0 ;1051 if (x < real::R_1()) 1052 return real::R_0(); 1003 1053 1004 1054 real ret = x; … … 1024 1074 *  if x == floor(x), return x 1025 1075 *  otherwise, return floor(x) + 1 */ 1026 if (x < real::R_0 )1076 if (x < real::R_0()) 1027 1077 return floor(x); 1028 1078 real ret = floor(x); … … 1030 1080 return ret; 1031 1081 else 1032 return ret + real::R_1 ;1082 return ret + real::R_1(); 1033 1083 } 1034 1084 1035 1085 template<> real round(real const &x) 1036 1086 { 1037 if (x < real::R_0 )1087 if (x < real::R_0()) 1038 1088 return round(x); 1039 1089 1040 return floor(x + (real::R_1 / 2));1090 return floor(x + (real::R_1() / 2)); 1041 1091 } 1042 1092 … … 1044 1094 { 1045 1095 if (!y) 1046 return real::R_0 ; /* FIXME: return NaN */1096 return real::R_0(); /* FIXME: return NaN */ 1047 1097 1048 1098 if (!x) … … 1057 1107 int switch_sign = x.m_signexp & 0x80000000u; 1058 1108 1059 real absx = fmod(fabs(x), real::R_PI * 2);1060 if (absx > real::R_PI )1061 { 1062 absx = real::R_PI ;1109 real absx = fmod(fabs(x), real::R_PI() * 2); 1110 if (absx > real::R_PI()) 1111 { 1112 absx = real::R_PI(); 1063 1113 switch_sign = !switch_sign; 1064 1114 } 1065 1115 1066 if (absx > real::R_PI_2 )1067 absx = real::R_PI  absx;1068 1069 real ret = real::R_0 , fact = real::R_1, xn = absx, mx2 = absx * absx;1116 if (absx > real::R_PI_2()) 1117 absx = real::R_PI()  absx; 1118 1119 real ret = real::R_0(), fact = real::R_1(), xn = absx, mx2 = absx * absx; 1070 1120 int i = 1; 1071 1121 for (;;) … … 1088 1138 template<> real cos(real const &x) 1089 1139 { 1090 return sin(real::R_PI_2  x);1140 return sin(real::R_PI_2()  x); 1091 1141 } 1092 1142 … … 1094 1144 { 1095 1145 /* Constrain input to [π,π] */ 1096 real y = fmod(x, real::R_PI );1146 real y = fmod(x, real::R_PI()); 1097 1147 1098 1148 /* Constrain input to [π/2,π/2] */ 1099 if (y < real::R_PI_2 )1100 y += real::R_PI ;1101 else if (y > real::R_PI_2 )1102 y = real::R_PI ;1149 if (y < real::R_PI_2()) 1150 y += real::R_PI(); 1151 else if (y > real::R_PI_2()) 1152 y = real::R_PI(); 1103 1153 1104 1154 /* In [π/4,π/4] return sin/cos */ 1105 if (fabs(y) <= real::R_PI_4 )1155 if (fabs(y) <= real::R_PI_4()) 1106 1156 return sin(y) / cos(y); 1107 1157 1108 1158 /* Otherwise, return cos/sin */ 1109 if (y > real::R_0 )1110 y = real::R_PI_2  y;1159 if (y > real::R_0()) 1160 y = real::R_PI_2()  y; 1111 1161 else 1112 y = real::R_PI_2  y;1162 y = real::R_PI_2()  y; 1113 1163 1114 1164 return cos(y) / sin(y); … … 1123 1173 * lose the precision around x=1. */ 1124 1174 real absx = fabs(x); 1125 int around_zero = (absx < (real::R_1 / 2));1175 int around_zero = (absx < (real::R_1() / 2)); 1126 1176 1127 1177 if (!around_zero) 1128 absx = sqrt((real::R_1  absx) / 2);1178 absx = sqrt((real::R_1()  absx) / 2); 1129 1179 1130 1180 real ret = absx, xn = absx, x2 = absx * absx, fact1 = 2, fact2 = 1; … … 1145 1195 1146 1196 if (around_zero) 1147 ret = is_asin ? ret : real::R_PI_2  ret;1197 ret = is_asin ? ret : real::R_PI_2()  ret; 1148 1198 else 1149 1199 { 1150 real adjust = is_negative ? real::R_PI : real::R_0;1200 real adjust = is_negative ? real::R_PI() : real::R_0(); 1151 1201 if (is_asin) 1152 ret = real::R_PI_2  adjust  ret * 2;1202 ret = real::R_PI_2()  adjust  ret * 2; 1153 1203 else 1154 1204 ret = adjust + ret * 2; … … 1194 1244 real absx = fabs(x); 1195 1245 1196 if (absx < (real::R_1 / 2))1246 if (absx < (real::R_1() / 2)) 1197 1247 { 1198 1248 real ret = x, xn = x, mx2 = x * x; … … 1210 1260 real ret = 0; 1211 1261 1212 if (absx < (real::R_3 / 2))1213 { 1214 real y = real::R_1  absx;1262 if (absx < (real::R_3() / 2)) 1263 { 1264 real y = real::R_1()  absx; 1215 1265 real yn = y, my2 = y * y; 1216 1266 for (int i = 0; ; i += 2) … … 1226 1276 yn *= my2; 1227 1277 } 1228 ret = real::R_PI_4  ret;1229 } 1230 else if (absx < real::R_2 )1231 { 1232 real y = (absx  real::R_SQRT3 ) / 2;1278 ret = real::R_PI_4()  ret; 1279 } 1280 else if (absx < real::R_2()) 1281 { 1282 real y = (absx  real::R_SQRT3()) / 2; 1233 1283 real yn = y, my2 = y * y; 1234 1284 for (int i = 1; ; i += 6) … … 1236 1286 real newret = ret + ((yn / (real)i) / 2); 1237 1287 yn *= y; 1238 newret = (real::R_SQRT3 / 2) * yn / (real)(i + 1);1288 newret = (real::R_SQRT3() / 2) * yn / (real)(i + 1); 1239 1289 yn *= y; 1240 1290 newret += yn / (real)(i + 2); 1241 1291 yn *= y; 1242 newret = (real::R_SQRT3 / 2) * yn / (real)(i + 3);1292 newret = (real::R_SQRT3() / 2) * yn / (real)(i + 3); 1243 1293 yn *= y; 1244 1294 newret += (yn / (real)(i + 4)) / 2; … … 1248 1298 yn *= my2; 1249 1299 } 1250 ret = real::R_PI_3 + ret;1300 ret = real::R_PI_3() + ret; 1251 1301 } 1252 1302 else … … 1263 1313 ret = newret; 1264 1314 } 1265 ret = real::R_PI_2  ret;1315 ret = real::R_PI_2()  ret; 1266 1316 } 1267 1317 … … 1278 1328 return y; 1279 1329 if (y.m_signexp >> 31) 1280 return real::R_PI ;1281 return real::R_PI ;1330 return real::R_PI(); 1331 return real::R_PI(); 1282 1332 } 1283 1333 … … 1285 1335 { 1286 1336 if (y.m_signexp >> 31) 1287 return real::R_PI ;1288 return real::R_PI ;1337 return real::R_PI(); 1338 return real::R_PI(); 1289 1339 } 1290 1340 … … 1292 1342 real z = y / x; 1293 1343 real ret = atan(z); 1294 if (x < real::R_0 )1295 ret += (y > real::R_0 ) ? real::R_PI : real::R_PI;1344 if (x < real::R_0()) 1345 ret += (y > real::R_0()) ? real::R_PI() : real::R_PI(); 1296 1346 return ret; 1297 1347 } … … 1334 1384 /* FIXME: better use int64_t when the cast is implemented */ 1335 1385 int exponent = ceil(log10(x)); 1336 x /= pow(R_10 , (real)exponent);1337 1338 if (x < R_1 )1339 { 1340 x *= R_10 ;1386 x /= pow(R_10(), (real)exponent); 1387 1388 if (x < R_1()) 1389 { 1390 x *= R_10(); 1341 1391 exponent; 1342 1392 } … … 1350 1400 *str++ = '.'; 1351 1401 x = real(digit); 1352 x *= R_10 ;1402 x *= R_10(); 1353 1403 } 1354 1404 … … 1379 1429 } 1380 1430 1381 template<> real const real::R_0 = (real)0.0;1382 template<> real const real::R_1 = (real)1.0;1383 template<> real const real::R_2 = (real)2.0;1384 template<> real const real::R_3 = (real)3.0;1385 template<> real const real::R_10 = (real)10.0;1386 1387 /*1388 * Initialisation order is important here:1389 *  fast_log() requires R_11390 *  log() requires R_LN21391 *  re() require R_21392 *  exp() requires R_0, R_1, R_LN21393 *  sqrt() requires R_31394 */1395 template<> real const real::R_LN2 = fast_log(R_2);1396 template<> real const real::R_LN10 = log(R_10);1397 template<> real const real::R_LOG2E = re(R_LN2);1398 template<> real const real::R_LOG10E = re(R_LN10);1399 template<> real const real::R_E = exp(R_1);1400 template<> real const real::R_PI = fast_pi();1401 template<> real const real::R_PI_2 = R_PI / 2;1402 template<> real const real::R_PI_3 = R_PI / R_3;1403 template<> real const real::R_PI_4 = R_PI / 4;1404 template<> real const real::R_1_PI = re(R_PI);1405 template<> real const real::R_2_PI = R_1_PI * 2;1406 template<> real const real::R_2_SQRTPI = re(sqrt(R_PI)) * 2;1407 template<> real const real::R_SQRT2 = sqrt(R_2);1408 template<> real const real::R_SQRT3 = sqrt(R_3);1409 template<> real const real::R_SQRT1_2 = R_SQRT2 / 2;1410 1411 1431 } /* namespace lol */ 1412 1432
Note: See TracChangeset
for help on using the changeset viewer.