llmath.h
上传用户:king477883
上传日期:2021-03-01
资源大小:9553k
文件大小:13k
源码类别:

游戏引擎

开发平台:

C++ Builder

  1. /** 
  2.  * @file llmath.h
  3.  * @brief Useful math constants and macros.
  4.  *
  5.  * $LicenseInfo:firstyear=2000&license=viewergpl$
  6.  * 
  7.  * Copyright (c) 2000-2010, Linden Research, Inc.
  8.  * 
  9.  * Second Life Viewer Source Code
  10.  * The source code in this file ("Source Code") is provided by Linden Lab
  11.  * to you under the terms of the GNU General Public License, version 2.0
  12.  * ("GPL"), unless you have obtained a separate licensing agreement
  13.  * ("Other License"), formally executed by you and Linden Lab.  Terms of
  14.  * the GPL can be found in doc/GPL-license.txt in this distribution, or
  15.  * online at http://secondlifegrid.net/programs/open_source/licensing/gplv2
  16.  * 
  17.  * There are special exceptions to the terms and conditions of the GPL as
  18.  * it is applied to this Source Code. View the full text of the exception
  19.  * in the file doc/FLOSS-exception.txt in this software distribution, or
  20.  * online at
  21.  * http://secondlifegrid.net/programs/open_source/licensing/flossexception
  22.  * 
  23.  * By copying, modifying or distributing this software, you acknowledge
  24.  * that you have read and understood your obligations described above,
  25.  * and agree to abide by those obligations.
  26.  * 
  27.  * ALL LINDEN LAB SOURCE CODE IS PROVIDED "AS IS." LINDEN LAB MAKES NO
  28.  * WARRANTIES, EXPRESS, IMPLIED OR OTHERWISE, REGARDING ITS ACCURACY,
  29.  * COMPLETENESS OR PERFORMANCE.
  30.  * $/LicenseInfo$
  31.  */
  32. #ifndef LLMATH_H
  33. #define LLMATH_H
  34. #include <cmath>
  35. #include <cstdlib>
  36. #include <complex>
  37. #include "lldefs.h"
  38. //#include "llstl.h" // *TODO: Remove when LLString is gone
  39. //#include "llstring.h" // *TODO: Remove when LLString is gone
  40. // lltut.h uses is_approx_equal_fraction(). This was moved to its own header
  41. // file in llcommon so we can use lltut.h for llcommon tests without making
  42. // llcommon depend on llmath.
  43. #include "is_approx_equal_fraction.h"
  44. // work around for Windows & older gcc non-standard function names.
  45. #if LL_WINDOWS
  46. #include <float.h>
  47. #define llisnan(val) _isnan(val)
  48. #define llfinite(val) _finite(val)
  49. #elif (LL_LINUX && __GNUC__ <= 2)
  50. #define llisnan(val) isnan(val)
  51. #define llfinite(val) isfinite(val)
  52. #elif LL_SOLARIS
  53. #define llisnan(val)    isnan(val)
  54. #define llfinite(val)   (val <= std::numeric_limits<double>::max())
  55. #else
  56. #define llisnan(val) std::isnan(val)
  57. #define llfinite(val) std::isfinite(val)
  58. #endif
  59. // Single Precision Floating Point Routines
  60. #ifndef fsqrtf
  61. #define fsqrtf(x) ((F32)sqrt((F64)(x)))
  62. #endif
  63. #ifndef sqrtf
  64. #define sqrtf(x) ((F32)sqrt((F64)(x)))
  65. #endif
  66. #ifndef cosf
  67. #define cosf(x) ((F32)cos((F64)(x)))
  68. #endif
  69. #ifndef sinf
  70. #define sinf(x) ((F32)sin((F64)(x)))
  71. #endif
  72. #ifndef tanf
  73. #define tanf(x) ((F32)tan((F64)(x)))
  74. #endif
  75. #ifndef acosf
  76. #define acosf(x) ((F32)acos((F64)(x)))
  77. #endif
  78. #ifndef powf
  79. #define powf(x,y) ((F32)pow((F64)(x),(F64)(y)))
  80. #endif
  81. const F32 GRAVITY = -9.8f;
  82. // mathematical constants
  83. const F32 F_PI = 3.1415926535897932384626433832795f;
  84. const F32 F_TWO_PI = 6.283185307179586476925286766559f;
  85. const F32 F_PI_BY_TWO = 1.5707963267948966192313216916398f;
  86. const F32 F_SQRT_TWO_PI = 2.506628274631000502415765284811f;
  87. const F32 F_E = 2.71828182845904523536f;
  88. const F32 F_SQRT2 = 1.4142135623730950488016887242097f;
  89. const F32 F_SQRT3 = 1.73205080756888288657986402541f;
  90. const F32 OO_SQRT2 = 0.7071067811865475244008443621049f;
  91. const F32 DEG_TO_RAD = 0.017453292519943295769236907684886f;
  92. const F32 RAD_TO_DEG = 57.295779513082320876798154814105f;
  93. const F32 F_APPROXIMATELY_ZERO = 0.00001f;
  94. const F32 F_LN2 = 0.69314718056f;
  95. const F32 OO_LN2 = 1.4426950408889634073599246810019f;
  96. const F32 F_ALMOST_ZERO = 0.0001f;
  97. const F32 F_ALMOST_ONE = 1.0f - F_ALMOST_ZERO;
  98. // BUG: Eliminate in favor of F_APPROXIMATELY_ZERO above?
  99. const F32 FP_MAG_THRESHOLD = 0.0000001f;
  100. // TODO: Replace with logic like is_approx_equal
  101. inline BOOL is_approx_zero( F32 f ) { return (-F_APPROXIMATELY_ZERO < f) && (f < F_APPROXIMATELY_ZERO); }
  102. // These functions work by interpreting sign+exp+mantissa as an unsigned
  103. // integer.
  104. // For example:
  105. // x = <sign>1 <exponent>00000010 <mantissa>00000000000000000000000
  106. // y = <sign>1 <exponent>00000001 <mantissa>11111111111111111111111
  107. //
  108. // interpreted as ints = 
  109. // x = 10000001000000000000000000000000
  110. // y = 10000000111111111111111111111111
  111. // which is clearly a different of 1 in the least significant bit
  112. // Values with the same exponent can be trivially shown to work.
  113. //
  114. // WARNING: Denormals of opposite sign do not work
  115. // x = <sign>1 <exponent>00000000 <mantissa>00000000000000000000001
  116. // y = <sign>0 <exponent>00000000 <mantissa>00000000000000000000001
  117. // Although these values differ by 2 in the LSB, the sign bit makes
  118. // the int comparison fail.
  119. //
  120. // WARNING: NaNs can compare equal
  121. // There is no special treatment of exceptional values like NaNs
  122. //
  123. // WARNING: Infinity is comparable with F32_MAX and negative 
  124. // infinity is comparable with F32_MIN
  125. inline BOOL is_approx_equal(F32 x, F32 y)
  126. {
  127. const S32 COMPARE_MANTISSA_UP_TO_BIT = 0x02;
  128. return (std::abs((S32) ((U32&)x - (U32&)y) ) < COMPARE_MANTISSA_UP_TO_BIT);
  129. }
  130. inline BOOL is_approx_equal(F64 x, F64 y)
  131. {
  132. const S64 COMPARE_MANTISSA_UP_TO_BIT = 0x02;
  133. return (std::abs((S32) ((U64&)x - (U64&)y) ) < COMPARE_MANTISSA_UP_TO_BIT);
  134. }
  135. inline S32 llabs(const S32 a)
  136. {
  137. return S32(std::labs(a));
  138. }
  139. inline F32 llabs(const F32 a)
  140. {
  141. return F32(std::fabs(a));
  142. }
  143. inline F64 llabs(const F64 a)
  144. {
  145. return F64(std::fabs(a));
  146. }
  147. inline S32 lltrunc( F32 f )
  148. {
  149. #if LL_WINDOWS && !defined( __INTEL_COMPILER )
  150. // Avoids changing the floating point control word.
  151. // Add or subtract 0.5 - epsilon and then round
  152. const static U32 zpfp[] = { 0xBEFFFFFF, 0x3EFFFFFF };
  153. S32 result;
  154. __asm {
  155. fld f
  156. mov eax, f
  157. shr eax, 29
  158. and eax, 4
  159. fadd dword ptr [zpfp + eax]
  160. fistp result
  161. }
  162. return result;
  163. #else
  164. return (S32)f;
  165. #endif
  166. }
  167. inline S32 lltrunc( F64 f )
  168. {
  169. return (S32)f;
  170. }
  171. inline S32 llfloor( F32 f )
  172. {
  173. #if LL_WINDOWS && !defined( __INTEL_COMPILER )
  174. // Avoids changing the floating point control word.
  175. // Accurate (unlike Stereopsis version) for all values between S32_MIN and S32_MAX and slightly faster than Stereopsis version.
  176. // Add -(0.5 - epsilon) and then round
  177. const U32 zpfp = 0xBEFFFFFF;
  178. S32 result;
  179. __asm { 
  180. fld f
  181. fadd dword ptr [zpfp]
  182. fistp result
  183. }
  184. return result;
  185. #else
  186. return (S32)floorf(f);
  187. #endif
  188. }
  189. inline S32 llceil( F32 f )
  190. {
  191. // This could probably be optimized, but this works.
  192. return (S32)ceil(f);
  193. }
  194. #ifndef BOGUS_ROUND
  195. // Use this round.  Does an arithmetic round (0.5 always rounds up)
  196. inline S32 llround(const F32 val)
  197. {
  198. return llfloor(val + 0.5f);
  199. }
  200. #else // BOGUS_ROUND
  201. // Old llround implementation - does banker's round (toward nearest even in the case of a 0.5.
  202. // Not using this because we don't have a consistent implementation on both platforms, use
  203. // llfloor(val + 0.5f), which is consistent on all platforms.
  204. inline S32 llround(const F32 val)
  205. {
  206. #if LL_WINDOWS
  207. // Note: assumes that the floating point control word is set to rounding mode (the default)
  208. S32 ret_val;
  209. _asm fld val
  210. _asm fistp ret_val;
  211. return ret_val;
  212. #elif LL_LINUX
  213. // Note: assumes that the floating point control word is set
  214. // to rounding mode (the default)
  215. S32 ret_val;
  216. __asm__ __volatile__( "flds %1    nt"
  217.   "fistpl %0  nt"
  218.   : "=m" (ret_val)
  219.   : "m" (val) );
  220. return ret_val;
  221. #else
  222. return llfloor(val + 0.5f);
  223. #endif
  224. }
  225. // A fast arithmentic round on intel, from Laurent de Soras http://ldesoras.free.fr
  226. inline int round_int(double x)
  227. {
  228. const float round_to_nearest = 0.5f;
  229. int i;
  230. __asm
  231. {
  232. fld x
  233. fadd st, st (0)
  234. fadd round_to_nearest
  235. fistp i
  236. sar i, 1
  237. }
  238. return (i);
  239. }
  240. #endif // BOGUS_ROUND
  241. inline F32 llround( F32 val, F32 nearest )
  242. {
  243. return F32(floor(val * (1.0f / nearest) + 0.5f)) * nearest;
  244. }
  245. inline F64 llround( F64 val, F64 nearest )
  246. {
  247. return F64(floor(val * (1.0 / nearest) + 0.5)) * nearest;
  248. }
  249. // these provide minimum peak error
  250. //
  251. // avg  error = -0.013049 
  252. // peak error = -31.4 dB
  253. // RMS  error = -28.1 dB
  254. const F32 FAST_MAG_ALPHA = 0.960433870103f;
  255. const F32 FAST_MAG_BETA = 0.397824734759f;
  256. // these provide minimum RMS error
  257. //
  258. // avg  error = 0.000003 
  259. // peak error = -32.6 dB
  260. // RMS  error = -25.7 dB
  261. //
  262. //const F32 FAST_MAG_ALPHA = 0.948059448969f;
  263. //const F32 FAST_MAG_BETA = 0.392699081699f;
  264. inline F32 fastMagnitude(F32 a, F32 b)
  265. a = (a > 0) ? a : -a;
  266. b = (b > 0) ? b : -b;
  267. return(FAST_MAG_ALPHA * llmax(a,b) + FAST_MAG_BETA * llmin(a,b));
  268. }
  269. ////////////////////
  270. //
  271. // Fast F32/S32 conversions
  272. //
  273. // Culled from www.stereopsis.com/FPU.html
  274. const F64 LL_DOUBLE_TO_FIX_MAGIC = 68719476736.0*1.5;     //2^36 * 1.5,  (52-_shiftamt=36) uses limited precisicion to floor
  275. const S32 LL_SHIFT_AMOUNT = 16;                    //16.16 fixed point representation,
  276. // Endian dependent code
  277. #ifdef LL_LITTLE_ENDIAN
  278. #define LL_EXP_INDEX 1
  279. #define LL_MAN_INDEX 0
  280. #else
  281. #define LL_EXP_INDEX 0
  282. #define LL_MAN_INDEX 1
  283. #endif
  284. /* Deprecated: use llround(), lltrunc(), or llfloor() instead
  285. // ================================================================================================
  286. // Real2Int
  287. // ================================================================================================
  288. inline S32 F64toS32(F64 val)
  289. {
  290. val = val + LL_DOUBLE_TO_FIX_MAGIC;
  291. return ((S32*)&val)[LL_MAN_INDEX] >> LL_SHIFT_AMOUNT; 
  292. }
  293. // ================================================================================================
  294. // Real2Int
  295. // ================================================================================================
  296. inline S32 F32toS32(F32 val)
  297. {
  298. return F64toS32 ((F64)val);
  299. }
  300. */
  301. ////////////////////////////////////////////////
  302. //
  303. // Fast exp and log
  304. //
  305. // Implementation of fast exp() approximation (from a paper by Nicol N. Schraudolph
  306. // http://www.inf.ethz.ch/~schraudo/pubs/exp.pdf
  307. static union
  308. {
  309. double d;
  310. struct
  311. {
  312. #ifdef LL_LITTLE_ENDIAN
  313. S32 j, i;
  314. #else
  315. S32 i, j;
  316. #endif
  317. } n;
  318. } LLECO; // not sure what the name means
  319. #define LL_EXP_A (1048576 * OO_LN2) // use 1512775 for integer
  320. #define LL_EXP_C (60801) // this value of C good for -4 < y < 4
  321. #define LL_FAST_EXP(y) (LLECO.n.i = llround(F32(LL_EXP_A*(y))) + (1072693248 - LL_EXP_C), LLECO.d)
  322. inline F32 llfastpow(const F32 x, const F32 y)
  323. {
  324. return (F32)(LL_FAST_EXP(y * log(x)));
  325. }
  326. inline F32 snap_to_sig_figs(F32 foo, S32 sig_figs)
  327. {
  328. // compute the power of ten
  329. F32 bar = 1.f;
  330. for (S32 i = 0; i < sig_figs; i++)
  331. {
  332. bar *= 10.f;
  333. }
  334. foo = (F32)llround(foo * bar);
  335. // shift back
  336. foo /= bar;
  337. return foo;
  338. }
  339. inline F32 lerp(F32 a, F32 b, F32 u) 
  340. {
  341. return a + ((b - a) * u);
  342. }
  343. inline F32 lerp2d(F32 x00, F32 x01, F32 x10, F32 x11, F32 u, F32 v)
  344. {
  345. F32 a = x00 + (x01-x00)*u;
  346. F32 b = x10 + (x11-x10)*u;
  347. F32 r = a + (b-a)*v;
  348. return r;
  349. }
  350. inline F32 ramp(F32 x, F32 a, F32 b)
  351. {
  352. return (a == b) ? 0.0f : ((a - x) / (a - b));
  353. }
  354. inline F32 rescale(F32 x, F32 x1, F32 x2, F32 y1, F32 y2)
  355. {
  356. return lerp(y1, y2, ramp(x, x1, x2));
  357. }
  358. inline F32 clamp_rescale(F32 x, F32 x1, F32 x2, F32 y1, F32 y2)
  359. {
  360. if (y1 < y2)
  361. {
  362. return llclamp(rescale(x,x1,x2,y1,y2),y1,y2);
  363. }
  364. else
  365. {
  366. return llclamp(rescale(x,x1,x2,y1,y2),y2,y1);
  367. }
  368. }
  369. inline F32 cubic_step( F32 x, F32 x0, F32 x1, F32 s0, F32 s1 )
  370. {
  371. if (x <= x0)
  372. return s0;
  373. if (x >= x1)
  374. return s1;
  375. F32 f = (x - x0) / (x1 - x0);
  376. return s0 + (s1 - s0) * (f * f) * (3.0f - 2.0f * f);
  377. }
  378. inline F32 cubic_step( F32 x )
  379. {
  380. x = llclampf(x);
  381. return (x * x) * (3.0f - 2.0f * x);
  382. }
  383. inline F32 quadratic_step( F32 x, F32 x0, F32 x1, F32 s0, F32 s1 )
  384. {
  385. if (x <= x0)
  386. return s0;
  387. if (x >= x1)
  388. return s1;
  389. F32 f = (x - x0) / (x1 - x0);
  390. F32 f_squared = f * f;
  391. return (s0 * (1.f - f_squared)) + ((s1 - s0) * f_squared);
  392. }
  393. inline F32 llsimple_angle(F32 angle)
  394. {
  395. while(angle <= -F_PI)
  396. angle += F_TWO_PI;
  397. while(angle >  F_PI)
  398. angle -= F_TWO_PI;
  399. return angle;
  400. }
  401. //SDK - Renamed this to get_lower_power_two, since this is what this actually does.
  402. inline U32 get_lower_power_two(U32 val, U32 max_power_two)
  403. {
  404. if(!max_power_two)
  405. {
  406. max_power_two = 1 << 31 ;
  407. }
  408. if(max_power_two & (max_power_two - 1))
  409. {
  410. return 0 ;
  411. }
  412. for(; val < max_power_two ; max_power_two >>= 1) ;
  413. return max_power_two ;
  414. }
  415. // calculate next highest power of two, limited by max_power_two
  416. // This is taken from a brilliant little code snipped on http://acius2.blogspot.com/2007/11/calculating-next-power-of-2.html
  417. // Basically we convert the binary to a solid string of 1's with the same
  418. // number of digits, then add one.  We subtract 1 initially to handle
  419. // the case where the number passed in is actually a power of two.
  420. // WARNING: this only works with 32 bit ints.
  421. inline U32 get_next_power_two(U32 val, U32 max_power_two)
  422. {
  423. if(!max_power_two)
  424. {
  425. max_power_two = 1 << 31 ;
  426. }
  427. if(val >= max_power_two)
  428. {
  429. return max_power_two;
  430. }
  431. val--;
  432. val = (val >> 1) | val;
  433. val = (val >> 2) | val;
  434. val = (val >> 4) | val;
  435. val = (val >> 8) | val;
  436. val = (val >> 16) | val;
  437. val++;
  438. return val;
  439. }
  440. //get the gaussian value given the linear distance from axis x and guassian value o
  441. inline F32 llgaussian(F32 x, F32 o)
  442. {
  443. return 1.f/(F_SQRT_TWO_PI*o)*powf(F_E, -(x*x)/(2*o*o));
  444. }
  445. #endif