TSMath.cs 17 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533
  1. using System;
  2. /* Copyright (C) <2009-2011> <Thorben Linneweber, Jitter Physics>
  3. *
  4. * This software is provided 'as-is', without any express or implied
  5. * warranty. In no event will the authors be held liable for any damages
  6. * arising from the use of this software.
  7. *
  8. * Permission is granted to anyone to use this software for any purpose,
  9. * including commercial applications, and to alter it and redistribute it
  10. * freely, subject to the following restrictions:
  11. *
  12. * 1. The origin of this software must not be misrepresented; you must not
  13. * claim that you wrote the original software. If you use this software
  14. * in a product, an acknowledgment in the product documentation would be
  15. * appreciated but is not required.
  16. * 2. Altered source versions must be plainly marked as such, and must not be
  17. * misrepresented as being the original software.
  18. * 3. This notice may not be removed or altered from any source distribution.
  19. */
  20. namespace TrueSync {
  21. /// <summary>
  22. /// Contains common math operations.
  23. /// </summary>
  24. public sealed class TSMath {
  25. /// <summary>
  26. /// PI constant.
  27. /// </summary>
  28. public static FP Pi = FP.Pi;
  29. /**
  30. * @brief PI over 2 constant.
  31. **/
  32. public static FP PiOver2 = FP.PiOver2;
  33. /// <summary>
  34. /// A small value often used to decide if numeric
  35. /// results are zero.
  36. /// </summary>
  37. public static FP Epsilon = FP.Epsilon;
  38. /**
  39. * @brief Degree to radians constant.
  40. **/
  41. public static FP Deg2Rad = FP.Deg2Rad;
  42. /**
  43. * @brief Radians to degree constant.
  44. **/
  45. public static FP Rad2Deg = FP.Rad2Deg;
  46. /**
  47. * @brief FP infinity.
  48. * */
  49. public static FP Infinity = FP.MaxValue;
  50. /// <summary>
  51. /// Gets the square root.
  52. /// </summary>
  53. /// <param name="number">The number to get the square root from.</param>
  54. /// <returns></returns>
  55. #region public static FP Sqrt(FP number)
  56. public static FP Sqrt(FP number) {
  57. return FP.Sqrt(number);
  58. }
  59. #endregion
  60. /// <summary>
  61. /// Gets the maximum number of two values.
  62. /// </summary>
  63. /// <param name="val1">The first value.</param>
  64. /// <param name="val2">The second value.</param>
  65. /// <returns>Returns the largest value.</returns>
  66. #region public static FP Max(FP val1, FP val2)
  67. public static FP Max(FP val1, FP val2) {
  68. return (val1 > val2) ? val1 : val2;
  69. }
  70. #endregion
  71. /// <summary>
  72. /// Gets the minimum number of two values.
  73. /// </summary>
  74. /// <param name="val1">The first value.</param>
  75. /// <param name="val2">The second value.</param>
  76. /// <returns>Returns the smallest value.</returns>
  77. #region public static FP Min(FP val1, FP val2)
  78. public static FP Min(FP val1, FP val2) {
  79. return (val1 < val2) ? val1 : val2;
  80. }
  81. #endregion
  82. /// <summary>
  83. /// Gets the maximum number of three values.
  84. /// </summary>
  85. /// <param name="val1">The first value.</param>
  86. /// <param name="val2">The second value.</param>
  87. /// <param name="val3">The third value.</param>
  88. /// <returns>Returns the largest value.</returns>
  89. #region public static FP Max(FP val1, FP val2,FP val3)
  90. public static FP Max(FP val1, FP val2, FP val3) {
  91. FP max12 = (val1 > val2) ? val1 : val2;
  92. return (max12 > val3) ? max12 : val3;
  93. }
  94. #endregion
  95. /// <summary>
  96. /// Returns a number which is within [min,max]
  97. /// </summary>
  98. /// <param name="value">The value to clamp.</param>
  99. /// <param name="min">The minimum value.</param>
  100. /// <param name="max">The maximum value.</param>
  101. /// <returns>The clamped value.</returns>
  102. #region public static FP Clamp(FP value, FP min, FP max)
  103. public static FP Clamp(FP value, FP min, FP max) {
  104. if (value < min)
  105. {
  106. value = min;
  107. return value;
  108. }
  109. if (value > max)
  110. {
  111. value = max;
  112. }
  113. return value;
  114. }
  115. #endregion
  116. /// <summary>
  117. /// Returns a number which is within [FP.Zero, FP.One]
  118. /// </summary>
  119. /// <param name="value">The value to clamp.</param>
  120. /// <returns>The clamped value.</returns>
  121. public static FP Clamp01(FP value)
  122. {
  123. if (value < FP.Zero)
  124. return FP.Zero;
  125. if (value > FP.One)
  126. return FP.One;
  127. return value;
  128. }
  129. /// <summary>
  130. /// Changes every sign of the matrix entry to '+'
  131. /// </summary>
  132. /// <param name="matrix">The matrix.</param>
  133. /// <param name="result">The absolute matrix.</param>
  134. #region public static void Absolute(ref JMatrix matrix,out JMatrix result)
  135. public static void Absolute(ref TSMatrix matrix, out TSMatrix result) {
  136. result.M11 = FP.Abs(matrix.M11);
  137. result.M12 = FP.Abs(matrix.M12);
  138. result.M13 = FP.Abs(matrix.M13);
  139. result.M21 = FP.Abs(matrix.M21);
  140. result.M22 = FP.Abs(matrix.M22);
  141. result.M23 = FP.Abs(matrix.M23);
  142. result.M31 = FP.Abs(matrix.M31);
  143. result.M32 = FP.Abs(matrix.M32);
  144. result.M33 = FP.Abs(matrix.M33);
  145. }
  146. #endregion
  147. /// <summary>
  148. /// Returns the sine of value.
  149. /// </summary>
  150. public static FP Sin(FP value) {
  151. return FP.Sin(value);
  152. }
  153. /// <summary>
  154. /// Returns the cosine of value.
  155. /// </summary>
  156. public static FP Cos(FP value) {
  157. return FP.Cos(value);
  158. }
  159. /// <summary>
  160. /// Returns the tan of value.
  161. /// </summary>
  162. public static FP Tan(FP value) {
  163. return FP.Tan(value);
  164. }
  165. /// <summary>
  166. /// Returns the arc sine of value.
  167. /// </summary>
  168. public static FP Asin(FP value) {
  169. return FP.Asin(value);
  170. }
  171. /// <summary>
  172. /// Returns the arc cosine of value.
  173. /// </summary>
  174. public static FP Acos(FP value) {
  175. return FP.Acos(value);
  176. }
  177. /// <summary>
  178. /// Returns the arc tan of value.
  179. /// </summary>
  180. public static FP Atan(FP value) {
  181. return FP.Atan(value);
  182. }
  183. /// <summary>
  184. /// Returns the arc tan of coordinates x-y.
  185. /// </summary>
  186. public static FP Atan2(FP y, FP x) {
  187. return FP.Atan2(y, x);
  188. }
  189. /// <summary>
  190. /// Returns the largest integer less than or equal to the specified number.
  191. /// </summary>
  192. public static FP Floor(FP value) {
  193. return FP.Floor(value);
  194. }
  195. /// <summary>
  196. /// Returns the smallest integral value that is greater than or equal to the specified number.
  197. /// </summary>
  198. public static FP Ceiling(FP value) {
  199. return value;
  200. }
  201. /// <summary>
  202. /// Rounds a value to the nearest integral value.
  203. /// If the value is halfway between an even and an uneven value, returns the even value.
  204. /// </summary>
  205. public static FP Round(FP value) {
  206. return FP.Round(value);
  207. }
  208. /// <summary>
  209. /// Returns a number indicating the sign of a Fix64 number.
  210. /// Returns 1 if the value is positive, 0 if is 0, and -1 if it is negative.
  211. /// </summary>
  212. public static int Sign(FP value) {
  213. return FP.Sign(value);
  214. }
  215. /// <summary>
  216. /// Returns the absolute value of a Fix64 number.
  217. /// Note: Abs(Fix64.MinValue) == Fix64.MaxValue.
  218. /// </summary>
  219. public static FP Abs(FP value) {
  220. return FP.Abs(value);
  221. }
  222. public static FP Barycentric(FP value1, FP value2, FP value3, FP amount1, FP amount2) {
  223. return value1 + (value2 - value1) * amount1 + (value3 - value1) * amount2;
  224. }
  225. public static FP CatmullRom(FP value1, FP value2, FP value3, FP value4, FP amount) {
  226. // Using formula from http://www.mvps.org/directx/articles/catmull/
  227. // Internally using FPs not to lose precission
  228. FP amountSquared = amount * amount;
  229. FP amountCubed = amountSquared * amount;
  230. return (FP)(0.5 * (2.0 * value2 +
  231. (value3 - value1) * amount +
  232. (2.0 * value1 - 5.0 * value2 + 4.0 * value3 - value4) * amountSquared +
  233. (3.0 * value2 - value1 - 3.0 * value3 + value4) * amountCubed));
  234. }
  235. public static FP Distance(FP value1, FP value2) {
  236. return FP.Abs(value1 - value2);
  237. }
  238. public static FP Hermite(FP value1, FP tangent1, FP value2, FP tangent2, FP amount) {
  239. // All transformed to FP not to lose precission
  240. // Otherwise, for high numbers of param:amount the result is NaN instead of Infinity
  241. FP v1 = value1, v2 = value2, t1 = tangent1, t2 = tangent2, s = amount, result;
  242. FP sCubed = s * s * s;
  243. FP sSquared = s * s;
  244. if (amount == 0f)
  245. result = value1;
  246. else if (amount == 1f)
  247. result = value2;
  248. else
  249. result = (2 * v1 - 2 * v2 + t2 + t1) * sCubed +
  250. (3 * v2 - 3 * v1 - 2 * t1 - t2) * sSquared +
  251. t1 * s +
  252. v1;
  253. return (FP)result;
  254. }
  255. public static FP Lerp(FP value1, FP value2, FP amount) {
  256. return value1 + (value2 - value1) * Clamp01(amount);
  257. }
  258. public static FP InverseLerp(FP value1, FP value2, FP amount) {
  259. if (value1 != value2)
  260. return Clamp01((amount - value1) / (value2 - value1));
  261. return FP.Zero;
  262. }
  263. public static FP SmoothStep(FP value1, FP value2, FP amount) {
  264. // It is expected that 0 < amount < 1
  265. // If amount < 0, return value1
  266. // If amount > 1, return value2
  267. FP result = Clamp(amount, 0f, 1f);
  268. result = Hermite(value1, 0f, value2, 0f, result);
  269. return result;
  270. }
  271. /// <summary>
  272. /// Returns 2 raised to the specified power.
  273. /// Provides at least 6 decimals of accuracy.
  274. /// </summary>
  275. internal static FP Pow2(FP x)
  276. {
  277. if (x.RawValue == 0)
  278. {
  279. return FP.One;
  280. }
  281. // Avoid negative arguments by exploiting that exp(-x) = 1/exp(x).
  282. bool neg = x.RawValue < 0;
  283. if (neg)
  284. {
  285. x = -x;
  286. }
  287. if (x == FP.One)
  288. {
  289. return neg ? FP.One / (FP)2 : (FP)2;
  290. }
  291. if (x >= FP.Log2Max)
  292. {
  293. return neg ? FP.One / FP.MaxValue : FP.MaxValue;
  294. }
  295. if (x <= FP.Log2Min)
  296. {
  297. return neg ? FP.MaxValue : FP.Zero;
  298. }
  299. /* The algorithm is based on the power series for exp(x):
  300. * http://en.wikipedia.org/wiki/Exponential_function#Formal_definition
  301. *
  302. * From term n, we get term n+1 by multiplying with x/n.
  303. * When the sum term drops to zero, we can stop summing.
  304. */
  305. int integerPart = (int)Floor(x);
  306. // Take fractional part of exponent
  307. x = FP.FromRaw(x.RawValue & 0x00000000FFFFFFFF);
  308. var result = FP.One;
  309. var term = FP.One;
  310. int i = 1;
  311. while (term.RawValue != 0)
  312. {
  313. term = FP.FastMul(FP.FastMul(x, term), FP.Ln2) / (FP)i;
  314. result += term;
  315. i++;
  316. }
  317. result = FP.FromRaw(result.RawValue << integerPart);
  318. if (neg)
  319. {
  320. result = FP.One / result;
  321. }
  322. return result;
  323. }
  324. /// <summary>
  325. /// Returns the base-2 logarithm of a specified number.
  326. /// Provides at least 9 decimals of accuracy.
  327. /// </summary>
  328. /// <exception cref="ArgumentOutOfRangeException">
  329. /// The argument was non-positive
  330. /// </exception>
  331. internal static FP Log2(FP x)
  332. {
  333. if (x.RawValue <= 0)
  334. {
  335. throw new ArgumentOutOfRangeException("Non-positive value passed to Ln", "x");
  336. }
  337. // This implementation is based on Clay. S. Turner's fast binary logarithm
  338. // algorithm (C. S. Turner, "A Fast Binary Logarithm Algorithm", IEEE Signal
  339. // Processing Mag., pp. 124,140, Sep. 2010.)
  340. long b = 1U << (FP.FRACTIONAL_PLACES - 1);
  341. long y = 0;
  342. long rawX = x.RawValue;
  343. while (rawX < FP.ONE)
  344. {
  345. rawX <<= 1;
  346. y -= FP.ONE;
  347. }
  348. while (rawX >= (FP.ONE << 1))
  349. {
  350. rawX >>= 1;
  351. y += FP.ONE;
  352. }
  353. var z = FP.FromRaw(rawX);
  354. for (int i = 0; i < FP.FRACTIONAL_PLACES; i++)
  355. {
  356. z = FP.FastMul(z, z);
  357. if (z.RawValue >= (FP.ONE << 1))
  358. {
  359. z = FP.FromRaw(z.RawValue >> 1);
  360. y += b;
  361. }
  362. b >>= 1;
  363. }
  364. return FP.FromRaw(y);
  365. }
  366. /// <summary>
  367. /// Returns the natural logarithm of a specified number.
  368. /// Provides at least 7 decimals of accuracy.
  369. /// </summary>
  370. /// <exception cref="ArgumentOutOfRangeException">
  371. /// The argument was non-positive
  372. /// </exception>
  373. public static FP Ln(FP x)
  374. {
  375. return FP.FastMul(Log2(x), FP.Ln2);
  376. }
  377. /// <summary>
  378. /// Returns a specified number raised to the specified power.
  379. /// Provides about 5 digits of accuracy for the result.
  380. /// </summary>
  381. /// <exception cref="DivideByZeroException">
  382. /// The base was zero, with a negative exponent
  383. /// </exception>
  384. /// <exception cref="ArgumentOutOfRangeException">
  385. /// The base was negative, with a non-zero exponent
  386. /// </exception>
  387. public static FP Pow(FP b, FP exp)
  388. {
  389. if (b == FP.One)
  390. {
  391. return FP.One;
  392. }
  393. if (exp.RawValue == 0)
  394. {
  395. return FP.One;
  396. }
  397. if (b.RawValue == 0)
  398. {
  399. if (exp.RawValue < 0)
  400. {
  401. //throw new DivideByZeroException();
  402. return FP.MaxValue;
  403. }
  404. return FP.Zero;
  405. }
  406. FP log2 = Log2(b);
  407. return Pow2(exp * log2);
  408. }
  409. public static FP MoveTowards(FP current, FP target, FP maxDelta)
  410. {
  411. if (Abs(target - current) <= maxDelta)
  412. return target;
  413. return (current + (Sign(target - current)) * maxDelta);
  414. }
  415. public static FP Repeat(FP t, FP length)
  416. {
  417. return (t - (Floor(t / length) * length));
  418. }
  419. public static FP DeltaAngle(FP current, FP target)
  420. {
  421. FP num = Repeat(target - current, (FP)360f);
  422. if (num > (FP)180f)
  423. {
  424. num -= (FP)360f;
  425. }
  426. return num;
  427. }
  428. public static FP MoveTowardsAngle(FP current, FP target, float maxDelta)
  429. {
  430. target = current + DeltaAngle(current, target);
  431. return MoveTowards(current, target, maxDelta);
  432. }
  433. public static FP SmoothDamp(FP current, FP target, ref FP currentVelocity, FP smoothTime, FP maxSpeed)
  434. {
  435. FP deltaTime = FP.EN2;
  436. return SmoothDamp(current, target, ref currentVelocity, smoothTime, maxSpeed, deltaTime);
  437. }
  438. public static FP SmoothDamp(FP current, FP target, ref FP currentVelocity, FP smoothTime)
  439. {
  440. FP deltaTime = FP.EN2;
  441. FP positiveInfinity = -FP.MaxValue;
  442. return SmoothDamp(current, target, ref currentVelocity, smoothTime, positiveInfinity, deltaTime);
  443. }
  444. public static FP SmoothDamp(FP current, FP target, ref FP currentVelocity, FP smoothTime, FP maxSpeed, FP deltaTime)
  445. {
  446. smoothTime = Max(FP.EN4, smoothTime);
  447. FP num = (FP)2f / smoothTime;
  448. FP num2 = num * deltaTime;
  449. FP num3 = FP.One / (((FP.One + num2) + (((FP)0.48f * num2) * num2)) + ((((FP)0.235f * num2) * num2) * num2));
  450. FP num4 = current - target;
  451. FP num5 = target;
  452. FP max = maxSpeed * smoothTime;
  453. num4 = Clamp(num4, -max, max);
  454. target = current - num4;
  455. FP num7 = (currentVelocity + (num * num4)) * deltaTime;
  456. currentVelocity = (currentVelocity - (num * num7)) * num3;
  457. FP num8 = target + ((num4 + num7) * num3);
  458. if (((num5 - current) > FP.Zero) == (num8 > num5))
  459. {
  460. num8 = num5;
  461. currentVelocity = (num8 - num5) / deltaTime;
  462. }
  463. return num8;
  464. }
  465. }
  466. }