InvSqrtEstimate.h 5.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160
  1. // David Eberly, Geometric Tools, Redmond WA 98052
  2. // Copyright (c) 1998-2020
  3. // Distributed under the Boost Software License, Version 1.0.
  4. // https://www.boost.org/LICENSE_1_0.txt
  5. // https://www.geometrictools.com/License/Boost/LICENSE_1_0.txt
  6. // Version: 4.0.2019.08.13
  7. #pragma once
  8. #include <Mathematics/Math.h>
  9. // Minimax polynomial approximations to 1/sqrt(x). The polynomial p(x) of
  10. // degree D minimizes the quantity maximum{|1/sqrt(x) - p(x)| : x in [1,2]}
  11. // over all polynomials of degree D.
  12. namespace WwiseGTE
  13. {
  14. template <typename Real>
  15. class InvSqrtEstimate
  16. {
  17. public:
  18. // The input constraint is x in [1,2]. For example,
  19. // float x; // in [1,2]
  20. // float result = InvSqrtEstimate<float>::Degree<3>(x);
  21. template <int D>
  22. inline static Real Degree(Real x)
  23. {
  24. Real t = x - (Real)1; // t in (0,1]
  25. return Evaluate(degree<D>(), t);
  26. }
  27. // The input constraint is x > 0. Range reduction is used to generate
  28. // a/ value y in [1,2], call Evaluate(y), and combine the output with
  29. // the proper exponent to obtain the approximation. For example,
  30. // float x; // x > 0
  31. // float result = InvSqrtEstimate<float>::DegreeRR<3>(x);
  32. template <int D>
  33. inline static Real DegreeRR(Real x)
  34. {
  35. Real adj, y;
  36. int p;
  37. Reduce(x, adj, y, p);
  38. Real poly = Degree<D>(y);
  39. Real result = Combine(adj, poly, p);
  40. return result;
  41. }
  42. private:
  43. // Metaprogramming and private implementation to allow specialization
  44. // of a template member function.
  45. template <int D> struct degree {};
  46. inline static Real Evaluate(degree<1>, Real t)
  47. {
  48. Real poly;
  49. poly = (Real)GTE_C_INVSQRT_DEG1_C1;
  50. poly = (Real)GTE_C_INVSQRT_DEG1_C0 + poly * t;
  51. return poly;
  52. }
  53. inline static Real Evaluate(degree<2>, Real t)
  54. {
  55. Real poly;
  56. poly = (Real)GTE_C_INVSQRT_DEG2_C2;
  57. poly = (Real)GTE_C_INVSQRT_DEG2_C1 + poly * t;
  58. poly = (Real)GTE_C_INVSQRT_DEG2_C0 + poly * t;
  59. return poly;
  60. }
  61. inline static Real Evaluate(degree<3>, Real t)
  62. {
  63. Real poly;
  64. poly = (Real)GTE_C_INVSQRT_DEG3_C3;
  65. poly = (Real)GTE_C_INVSQRT_DEG3_C2 + poly * t;
  66. poly = (Real)GTE_C_INVSQRT_DEG3_C1 + poly * t;
  67. poly = (Real)GTE_C_INVSQRT_DEG3_C0 + poly * t;
  68. return poly;
  69. }
  70. inline static Real Evaluate(degree<4>, Real t)
  71. {
  72. Real poly;
  73. poly = (Real)GTE_C_INVSQRT_DEG4_C4;
  74. poly = (Real)GTE_C_INVSQRT_DEG4_C3 + poly * t;
  75. poly = (Real)GTE_C_INVSQRT_DEG4_C2 + poly * t;
  76. poly = (Real)GTE_C_INVSQRT_DEG4_C1 + poly * t;
  77. poly = (Real)GTE_C_INVSQRT_DEG4_C0 + poly * t;
  78. return poly;
  79. }
  80. inline static Real Evaluate(degree<5>, Real t)
  81. {
  82. Real poly;
  83. poly = (Real)GTE_C_INVSQRT_DEG5_C5;
  84. poly = (Real)GTE_C_INVSQRT_DEG5_C4 + poly * t;
  85. poly = (Real)GTE_C_INVSQRT_DEG5_C3 + poly * t;
  86. poly = (Real)GTE_C_INVSQRT_DEG5_C2 + poly * t;
  87. poly = (Real)GTE_C_INVSQRT_DEG5_C1 + poly * t;
  88. poly = (Real)GTE_C_INVSQRT_DEG5_C0 + poly * t;
  89. return poly;
  90. }
  91. inline static Real Evaluate(degree<6>, Real t)
  92. {
  93. Real poly;
  94. poly = (Real)GTE_C_INVSQRT_DEG6_C6;
  95. poly = (Real)GTE_C_INVSQRT_DEG6_C5 + poly * t;
  96. poly = (Real)GTE_C_INVSQRT_DEG6_C4 + poly * t;
  97. poly = (Real)GTE_C_INVSQRT_DEG6_C3 + poly * t;
  98. poly = (Real)GTE_C_INVSQRT_DEG6_C2 + poly * t;
  99. poly = (Real)GTE_C_INVSQRT_DEG6_C1 + poly * t;
  100. poly = (Real)GTE_C_INVSQRT_DEG6_C0 + poly * t;
  101. return poly;
  102. }
  103. inline static Real Evaluate(degree<7>, Real t)
  104. {
  105. Real poly;
  106. poly = (Real)GTE_C_INVSQRT_DEG7_C7;
  107. poly = (Real)GTE_C_INVSQRT_DEG7_C6 + poly * t;
  108. poly = (Real)GTE_C_INVSQRT_DEG7_C5 + poly * t;
  109. poly = (Real)GTE_C_INVSQRT_DEG7_C4 + poly * t;
  110. poly = (Real)GTE_C_INVSQRT_DEG7_C3 + poly * t;
  111. poly = (Real)GTE_C_INVSQRT_DEG7_C2 + poly * t;
  112. poly = (Real)GTE_C_INVSQRT_DEG7_C1 + poly * t;
  113. poly = (Real)GTE_C_INVSQRT_DEG7_C0 + poly * t;
  114. return poly;
  115. }
  116. inline static Real Evaluate(degree<8>, Real t)
  117. {
  118. Real poly;
  119. poly = (Real)GTE_C_INVSQRT_DEG8_C8;
  120. poly = (Real)GTE_C_INVSQRT_DEG8_C7 + poly * t;
  121. poly = (Real)GTE_C_INVSQRT_DEG8_C6 + poly * t;
  122. poly = (Real)GTE_C_INVSQRT_DEG8_C5 + poly * t;
  123. poly = (Real)GTE_C_INVSQRT_DEG8_C4 + poly * t;
  124. poly = (Real)GTE_C_INVSQRT_DEG8_C3 + poly * t;
  125. poly = (Real)GTE_C_INVSQRT_DEG8_C2 + poly * t;
  126. poly = (Real)GTE_C_INVSQRT_DEG8_C1 + poly * t;
  127. poly = (Real)GTE_C_INVSQRT_DEG8_C0 + poly * t;
  128. return poly;
  129. }
  130. // Support for range reduction.
  131. inline static void Reduce(Real x, Real& adj, Real& y, int& p)
  132. {
  133. y = std::frexp(x, &p); // y in [1/2,1)
  134. y = ((Real)2) * y; // y in [1,2)
  135. --p;
  136. adj = (1 & p) * (Real)GTE_C_INV_SQRT_2 + (1 & ~p) * (Real)1;
  137. p = -(p >> 1);
  138. }
  139. inline static Real Combine(Real adj, Real y, int p)
  140. {
  141. return adj * std::ldexp(y, p);
  142. }
  143. };
  144. }